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1.  INTRODUCTION  AND  RESEARCH  OVERVIEW 


1.1  Dynamic  Hydrology  and  Complex  Terrain 

The  focus  of  this  research  has  been  the  study  of  spatio-temporal  dynamics  of  the 
rainfall-runoff  process  in  the  common  situation  where  topographic  relief  and  shallow 
subsurface  moisture  control  the  storage  and  timing  of  runoff  on  the  landscape.  The 
geometry  (shape,  curvature,  etc.)  of  the  hillslope,  catchment  or  river  basin  exerts  a  major 
control  on  the  hydrologic  storage-response,  by  defining  the  domain  and  boundary 
conditions  of  moisture  storage.  The  difficulties  stem  from  the  fact  that  the  governing 
equation  known  as  Richards  equation  (Richards,  1931)  is  strongly  nonlinear,  requiring  the 
solution  of  extremely  large  systems  of  equations  even  for  small  problems,  and  observed 
data  (eg.  soil  moisture,  groundwater  levels,  soil  and  geologic  permeability)  are  typically  not 
well  resolved  in  space  and/or  time.  Although  the  form  of  the  nonlinearities  in  the 
governing  local  equations  may  be  mathematically  understood,  the  expression  of  these 
nonlinearities  in  the  integrated  spatio-temporal  response  of  hillslope  and  watershed  runoff, 
are  not  understood  at  present  Topographic,  soil  and  geologic  spatial  variation  further  adds 
to  the  computational  burden  of  solving  the  four  dimensional  Richards  equation.  Modeling 
the  physical  processes  of  such  a  system  over  multiple  spatial  scales  remains  a  fundamental 
challenge  for  hydrologic  science. 

The  role  of  topographic  structure  on  hydrologic  processes  and  modeling  is  a 
fundamental  theme  of  this  research  effort.  The  research  will  attempt  to  incorporate  the 
essentials  of  this  structure  directly  into  the  model  fonnulation.  This  approach  is  felt  to  be 
particularly  useful  where  hydrologic  predictions  must  be  made  over  multiple  scales 
(theories)  appropriate  to  hillslopes,  watersheds  and  river  basins.  In  this  regard  we  make 
the  following  assumptions:  (1)  The  hillslope  is  the  fundamental  spatial  scale  characterizing 
a  storage  reservior  which  produces  runoff.  The  hillslope  itself  is  defined  as  that  region 
bounded  above  by  a  topographic/hydrologic  divide,  bounded  laterally  by  trajectories 
following  the  steepest  path  from  divide  to  stream,  and  bounded  below  by  a  line  or  surface 
of  discharge.  It  is  noted  that  the  boundary  conditions  of  all  hillslopes  are  directly  tied  to  the 
drainage  network.  (2)  The  classical  hypsometric,  or  area-elevation  relationship,  represents 
the  physical  domain  of  a  hydrologic  response  unit.  As  such,  it  can  be  applied  to  all  scales 
of  interest  (e.g.  the  hillslope,  watershed  or  river  basin  scale).  (3)  The  hypsometric  relation 
F((|))  is  defined  as  a  probability  distribution  in  the  following  sense.  F(<t>)  is  the  fraction  of 
area  below  an  elevation  contour  (|),  normalized  by  the  total  area.  Likewise  there  exists  a 
complementary  hypsometric  relationship  F(\i/),  where  V|f  represent  the  orthogonal 
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trajectories  of  <|i.  F(t(/)  is  defined  as  the  fraction  of  area  contained  below  a  y  contour.  (4) 
Spatial  averaging  of  hydrologic  state  variables  (e.g.  soil  misture,  saturated  groundwater 
levels,  etc)  should  be  made  with  respect  to  the  particular  hydrologic  unit  The  appropriate 
weighting  function  for  spatial  averaging  is  the  hypsometric  distribution.  Chapter  2 
discusses  in  detail  the  implications  of  hypsometry  to  modeling  and  data  structures  within  a 
hydrologic  response  unit. 

A  major  effort  of  the  last  three  years  has  been  to  devise  a  methodology  for 
constructing  spatially  integrated  constitutive  relations  for  storage  and  flux  of  water  on 
hillslopes  and  small  watersheds.  As  the  field  data  is  typically  sparse,  the  analysis  must  be 
supported  by  numerical  experiments.  A  comprehensive  series  of  numerical  experiments 
were  conducted  as  part  of  this  research,  and  spatially  integrated  constitutive  relations 
required  in  the  dynamical  model  were  proposed.  In  Chapter  3,  steady-state  flux-storage 
relationships  are  constructed  from  Richards  equation  subject  to  variations  in  slope  shape, 
slope  gradient,  aspect  ratio,  and  soil  type.  Steady-state  flux-storage  relationships  are 
parameterized  using  a  polynomial  relationship.  Following  the  steady-state  analysis,  the 
dynamic  storage-flux  relationships  are  examined  by  a  series  of  numerical  experiments  with 
a  range  of  input  forcing  pattern,  slope  shape,  soil  type,  and  boundary  conditions. 

In  chapter  4  the  research  will  explore  the  storage-flux  relationships  for  Watershed 
No.  2  of  the  Shale  Hills  Watersheds  in  central  PA  using  field  measurements  of  soil 
moisture,  water  table  elevation,  and  streamflow  (Lynch,  1976),  It  is  shown  that  the 
hypsometric  duistribution  of  the  hillslope,  and  spatial  scaling  transformations  of  the  surface 
and  bedrock  topography,  allow  construction  of  integrated  state  variables  of  the  hillslope. 
Namely  the  soil  moisture  y  and  saturated  storage  q.  Storage-flux  relationships  are  then 
established  for  the  watershed.  The  analysis  shows  that  the  dynamical  behaviour  between  y 
and  q  exhibits  an  inverse  relationship,  while  storage-flux  relationships  display  hysteretic 
behaviour.  The  field  analysis  demonstrates  that  local  topographic  scaling  is  a  useful  tool 
for  modeling  where  sparse  and  scattered  field  data  exists.  Chapter  5  discusses  solutions  to 
the  two-state  dynamical  system,  and  Chapter  6  examines  scaling  relations  of  digital 
elevation  models,  and  the  application  of  mapping  seasonal  precipitation  and  evaporation  in 
mountainous  terrain. 

To  summarize,  the  research  uses  a  hypsometric  basis  to  construct  spatial  averages 
of  hydrologic  state  variables  at  the  hillslope  and  watershed  scale.  The  basis  is  applied  to 
the  partial  differential  equations  for  unsaturated-saturated  flow  to  arrive  at  a  low¬ 
dimensional  predictive  model  for  rainfall-storage-runoff,  and  to  the  interpretation  of  field 
data  at  Shale  Hills  watershed  No.  2  to  test  the  concept.  Numerical  experiments  are  used  to 
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define  the  form  of  integrated  constitutive  relations  as  a  function  of  hillsope  geraietiy  (relief, 
slope,  curvature,  convergence,  soil  type,  soil  depth). 

1.2  A  Brief  Review  of  the  Literature 

Beven  (1983)  provides  a  comprehensive  review  of  the  role  of  basin  structure  on 
lunoff  generation.  He  points  out  the  complicating  issues  of  spatial  and  temporal  variability. 
Nonetheless,  he  contends  that  the  'ordered  nature'  of  drainage  basins  suggests  it  should  be 
possible  to  develop  a  unifying  theory  relating  basin  structure  to  runoff  response.  Siviplan  et 
al.  (1990)  provide  an  excellent  overview  of  the  importance  of  hydrologic  similarity  and 
scale  on  the  problem  of  rainfall-runoff.  The  concept  of  a  geomorphic  unit  hydrograph,  of 
Rodriguez-Iturbe  and  Valdes  (1979),  and  Gupta  et  al.  (1980)  provides  a  fundamental  link 
between  the  stream  network,  stream  order,  contributing  area,  and  empirical  models  of 
hydrologic  response.  The  present  research  serves  to  extend  the  concept  of  a 
hydrogeomorphic  response  unit  to  include  a  physically-based  and  nonlinear  response  at  the 
hillslope  scale. 

On  the  experimental  side,  Dunne  (1978)  summarizes  the  field  evidence  for  currently 
accepted  physical  mechanisms  generating  streamflow.  He  suggests  that  the  various 
processes  (Horton  overland  flow,  subsurface  flow,  and  saturation  overland  flow)  are  not 
mutually  exclusive,  but  rather  are  synergistic,  where  multiple  pathways  of  storm  runoff 
depend  on  physical  geography,  initial  condition,  and  intensity  of  the  hydrologic  event. 
Jeppson  and  Shreiber  (1972),  using  early  data  from  the  Upper  Sheep  Creek  experiment, 
were  one  of  the  first  to  apply  two-dimensional  unsaturated  flow  theory  to  hillslope 
modeling.  Subsequently  Stephenson  and  Freeze  (1974)  simulated  the  Upper  Sheep 
hillslope  as  a  two-dimensional,  transient  system  using  finite  differences.  During  the  70's, 
many  advances  were  made  in  algorithms  for  numerical  computation  in  subsurface 
hydrology  (Freeze,  1971;  Neuman,  Feddes  and  Bresler  1975),  and  in  instrumentation  and 
experimentation  of  soil  water  processes  in  the  field.  Conceptual  understanding  of  soil- 
water  processes  on  natural  landscapes  such  as  macropore  flow,  spatial  variability,  basin 
structure  etc.,  was  a  major  focus  of  research  during  this  period  (see  Beven,  1983,  for 
review).  At  the  basin  scale,  Klemes  (1983)  examined  the  monthly  dynamic  relationship 
between  runoff  and  basin  storage  using  rainfall  and  runoff  data  from  the  Rainy  Lake  Basin 
in  western  Ontario,  Canada.  The  storage-runoff  relationship  showed  a  complex  pattern, 
however  a  linear  relationship  was  used  based  on  the  winter  season  data.  In  a  small 
hillslope  catchment  near  the  western  suburbs  of  Tokyo,  Ando  et  al.  (1983)  measured 
groundwater  level  and  streamflow.  Using  these  data,  a  quadratic  relationship  between 


3 


groundwater  runoff  and  groundwater  storage  was  derived.  On  the  basis  of  this 
relationship,  hourly  and  daily  time  scale  hydrologic  models  were  developed. 

Recently  the  dynamics  of  soil  moisture  over  large  regions  has  begun  to  receive 
important  attention.  Rodriguez-Iturbe  et  ai.  (1991)  have  mined  the  nonlinear  dynamics 
of  soil  moisture  at,  what  they  refer  to  as  the  climate-scale.  In  this  research  a  statistical- 
dynamical  model  is  developed  to  represent  surface  hydrology  over  continental  regions. 
Famiglietti  and  Wood  (1991)  specifically  deal  with  development  of  a  catchment-scale  water 
balance  model  for  an  intensively  monitored  field  experiment  in  Kansas. 

Klemes  (1983)  emphasizes  testability  and  observability  as  important  elements  of 
conceptualization  and  model  development.  He  argues  that  many  conceptual  hydrologic 
models  are  developed  without  considering  their  testability.  Stephenson  and  Freeze  (1974) 
point  out  the  difficulties  of  calibrating  and  testing  a  physically  based  model  (i.e.,  Richards 
equation),  even  at  a  well-equipped  experimental  watershed.  The  problem  of  what  to 
measure,  and  where  to  measure  it  remains  an  important  issue.  The  application  of 
dynamical  systems  methodology  as  an  organizing  framework  to  study  hydrologic 
processes  in  complex  terrain,  offers  a  consistent  way  of  linking  data  to  models,  and 
ultimately  to  satisfy  requirements  of  observability  and  testiblility. 
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*•.  INTEGRAL  BALANCE  MODELS  AND  TOPOGRAPHIC  STRUCTURE 


2.1  The  Integral-Balance  and  State-Space  Concent 

It  can  be  said  that  the  fundamental  equation  of  hydrology  is  represented  by  the 
simple  water  balance  or  conservation  of  fluid  volume 

^  =  I-Q  (1) 

dt 

where  S  represents  the  fluid  storage  volume,  I  the  volumetric  input  and  Q  the  volumetric 
discharge  or  output.  For  a  hydrologic  system  with  interacting  physical  processes  and  time 
scales,  (1)  is  represented  in  the  general  state-space  form 

^  =  u(it)-/(tt)  (2) 

dt 

where  storage  is  now  a  vector  of  state  variables  %  =  input  flux  vector 

u  =  {u,,U2,...u„}  ,  and  the  output  flux  vector  /  =  {/p/j—  Zn)-  Formulating  the  particular 
dynamical  model,  requires  a  determination  of  the  geometry  and  spatial  scale  of  the  system, 
the  corresponding  dimension  of  the  state-space  n  and  finally,  estimation  of  the  particular 
form  of  the  constitutive  or  storage-flux  relationship  /(^,t),  u(|,t)  at  the  particular  scale. 
Clearly,  as  scales  and  geometry  change,  the  dimension  of  the  state-space  may  change  and 
the  storage-flux  relationships  may  change  as  well.  Nonetheless,  the  state-space  approach 
offers  a  general  framework  to  study  the  dynamics  of  hydrologic  systems  across  multiple 
scales. 

Throughout  this  research  it  is  assumed  that  state  variables  (e.g.  soil  moisture  and 
saturated  storage)  can  be  represented  by  spatial  averages.  These  averages  must  taken  with 
respect  to  the  geometry  of  the  hydrologic  flow  system  in  question,  and  not  simply  averaged 
over  an  arbitrary  region.  In  this  section  it  is  shown  that  the  hypsometric  or  area-elevation 
relation,  when  redefined  in  terms  of  probability,  serves  as  a  physical  meaningful  weighting 
function. 

To  proceed,  we  consider  a  general  space-time  representation  for  the  state  variable 
y(x,t)  of  the  form 

m 

y(x,t)  =  ^yk(t)ak(x)  (3) 

k-l 

where  yk(t)  represents  the  kth  amplitude  and  a^fx)  are  a  corresponding  set  of  spatial 
modes  or  shape  functions  with  respect  to  the  position  vector  x={xi,X2,X3}.  The  next 
section  develops  the  relationship  between  the  a^(x),  hypsometry  and  spatial  averages. 
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2.2  Hypsometry,  Probability  and  Topographic  Shape  Functions 

A  classical  descriptor  of  the  morphology  of  a  watershed  or  river  basin  uses  the 
concept  of  the  hypsometric  function  or  area-elevation  relationship,  introduced  to  the 
hyrology  literature  by  Langbein  et  al  (1947).  The  hypsometric  curve  is  an  empirical, 
nondimensional  curve  relating  the  fraction  of  catchment  area  contained  by  an  elevation 
contour.  It  is  defined  for  any  scale  of  hydrologic  response,  from  hillslope,  to  watershed, 
or  river  basin.  The  hypsometric  curve  or  its  integral  has  been  shown  to  be  related  to 
geomorphic  characteristics  of  the  drainage  network.  For  example  Hack  (1957)  relates 
drainage  basin  area  to  the  length  of  longest  channel.  Geologic  interpretations  of  the 
characteristic  shapes  of  the  hypsometric  curve  for  a  given  landform  (see  Schumm  (1956), 
or  Sheideggar  (1970)),  generally  assert  that  the  present  state  of  the  landform  depends  on 
the  balance  of  tectonic  and  erosive  activity,  and  degree  of  disequilibrium.  Notwithstanding 
these  implications,  the  hypsometric  relation  is  fundamentally  a  means  to  reduce  the  spatial 
dimension  of  a  topographic  surface  to  a  simple  function  of  elevation.  In  the  context  of 
hydrology,  the  hypsometric  relation  defines  the  domain  over  which  state  variables  (e.g.  soil 
moisture  and  saturated  storage)  contribute  to  runoff. 

The  hypsometry  of  a  hydrologic  unit  is  defined 

^  =  1-F«1>)  (4) 

where  ((>=z  represents  the  normalized  elevation  contour  above  base  level  (0<z^l),  A(<1>)  is 
area  associated  with  zx|>,  Aj  is  the  total  area  of  the  basin,  and  F(<|>)  is  the  fraction  of  total 
area  enclosed  by  z<((>  in  the  catchment,  or  in  terms  of  probability,  F(()))=prob(z<()>). 
Likewise,  we  may  define  a  function  orthogonal  to  (t>,  describing  the  path  or  trajectory  of  the 
landform,  \|r.  Both  0  and  y  are  defined  with  respect  to  the  natural  coordinates  of  the 
landform  {s,n}.  Figure  2.1  illustrates  the  elevation  contours  and  orthogonal  trajectories  for 
a  small  watershed  near  Penn  State  campus  known  as  Shale  Hills  Experimental  Watershed 
(J.  Lynch,  1976).  Also  indicated  in  Figure  1  is  the  hypsometric  relation  defined  by  (4). 
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Figure  2.1  Elevation  contours  and  orthogonal  trajectories  ^  and  y,  in  the  natural  coordinate 
system  {s,  n)  for  the  Shale  Hills  Watershed  #2.  Also  shown  is  the  hypsometric  relation. 


The  rationale  for  a  probabilistic  interpretation  of  the  hypsometric  relation,  is  to 
construct  a  weighting  function  for  spatial  averaging  of  state  variables.  For  the  state  variable 
y(<|>),  the  kth  spatial  moment  is  given 

EIy^(|»)l=  jy^iWffW;  .  (5). 

t«0 


Practical  experience  suggests  that  local  effects  of  aspect,  vegetation,  soils,  geology, 
etc.  may  influence  the  y-process.  In  these  cases  it  will  be  useful  to  construct  conditional 
spatial  averages  over  subregions  within  the  watershed.  A  subregion  is  defined  as  the  area 
contained  by  a  pair  of  flow  trajectories,  A\l^=\^f2*Vl  which  originate  at  the  local  topographic 
divide,  and  intersect  the  stream  at  local  base  level,  llie  conditional  distribution  for  the 
subregion  with  elevation  change  zi«t><Z2  is  given  by 


F((|llAv)= 

J  F(Z2)-F(z,) 

=  0.  <l><z, 

=  1,  4>  >  Zj  ’ 


and  the  conditional  density  is  defined 


d(|> 


ff(<j))d(l)  E(z2)“F(Zi) 


J*! 


(6) 


(7) 


Conditional  expectation  of  the  state  variable,  y(<{>)  within  the  subregion  Ay  is  given  by 
E[y((j»)IAv]  =  f  y(<l))  f(<|)l  Ay)d<{)  ( 8 ) 


If  the  catchment  is  made  up  of  non-overlapping  subregions  i=l,2,...m  the  density 
for  the  overall  region  is 


(9) 


1  a 


where  fi(0)  is  now  the  local  density  and  Xi  the  fraction  of  area  associated  with  Ayj.  The 
^atial  average  for  the  subregion  is  given  by 


Yi  =  Elyi(0)iAv 


«i))f.«l>lAv.)d(> 


The  spatial  average  for  the  entire  catchment  is  then  determined  from  (10)  and 


m 


(10) 


(11) 


A  general  approach  to  the  problem  of  hypsometric  subregions,  can  be  devised  by 
looking  at  the  joint  distribution  of  elevation  <{»  and  the  orthogonal  trajectories  y.  Let  Fy(vr) 
be  the  fraction  of  catchment  area  below  a  given  V|r. 

F^(¥)=  (12) 

total  area 

FK<1>)  is  now  given  in  complementary  form 

c  _ ^A(z)^^  /  1  \ 


(12) 


F*(0)  = 


total  area 


(13) 


The  relative  frequency  or  specific  area  for  each  is  given  by 

dF.(6)  dF.  -v) 

f  f  =  (14) 

^  dy 

The  functions  <j)  and  \|r  are  orthogonal,  and  if  they  also  happen  to  be  statistically 
independent ,  the  joint  distribution  can  be  written 

=  ^(V)  (15) 

Applying  the  usual  rules  for  determining  the  marginal  density  from  the  joint  density  yields 


^(<l>)  =  jF#»(<l>.V)dv;  fy(¥)  =  JF♦v(<^.V)d<t>; 


The  conditional  density  for  <(>  is 


(16) 


=  =  (17) 

To  estimate  spatial  moments  of  the  process  y=g(<i>,¥)>  we  take  the  joint-conditional 
expectation  over  the  subregion  A(t>  and  A¥ 

Vj*j 

E[yMA0,Av]=  JJ  g'‘(4>.V)fy,^(«l>IA<)»  fy,Ay(¥lAv)d(t>d¥.  ( 18) 
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Figure  2.2.  illustrates  the  complementary  functions  and  Fy(>|r)  for  the  Shale  Hills 
watershed  of  Figure  2.1,  and  indicates  a  subregion  over  which  spatial  moments  are  locally 
conditioned  by  discrete  pairs  of  elevation  contours  and  flow  paths. 

1.0 
0.8 
0.6 
0.4 
0.2 
0.0 

0.0  04  0.8 

Elevation 

Figure  2.2  The  complementary  hypsometric  distributions  F^((j)),  the  fraction  area  enclosed 
by  and  Fy(\ir),  the  area  enclosed  by  \j;<z.  The  hillslope  (shaded  region)  represents  the 
conditional  area  fraction  bounded  by  the  region  A0  and  Ay. 

The  above  probability  arguments  are  easily  extended  to  the  general  the  spatio- 
temporal  process  represented  by  (1).  For  the  special  case  where  a  hypsometric  basis  is 
appropriate  ( x-  >  4) .  It  is  instructive  to  decompose  y((j),t)  into  a  time  average  ( y)  and  a 
stationary  time-fluctuation  (y') 

y(<t>.t)=:y((|>)+y'(<|>,t)  (19) 

where  y(<|»)  is  a  simple  function  of  elevation  (<})),  while  the  perturbation  y'  ((j),t)  varies  with 
elevation  and  time.  Now  let  y'(<l>it),  have  the  expansion 

y'(<|).t)  =  X  yk(0a,(<l))  (20) 

k 

where  yk(t)  is  a  zero  mean  stationary  process  in  time  representing  stationary  harmonic 
and/or  random  contributions  to  the  overall  process,  and  tXk((}>)  represents  the  basis 
functions  or  shape  functions  for  the  same  component  of  the  y  fluctuation.  The  y-process 
now  has  the  representation 

y(<l».t)  =  y(<|))+Xyk(Oak(<{>)  (21) 

k 

Taking  the  first  moment  with  respect  to  elevation  ({>  yields 

M(„(t)  =  EJy((|),t)]  =  EJy(<)))l-»-y,(t)E,[a,((l»]4-yj(t)E[a2((l))]-»-.... 

=  J  y(<l>)f*(<l»)d<f>+  yi(t)J  a,(<|>)f,(<|>)d(|>+ yj(t)J  aj«|>)f,((|>)d(|n-. 

♦  ♦  ♦ 

=  y+Piyi(0+P2y2(0+. 


(22) 


The  coefficients  Pk  represents  the  elevation-weighted  average  contribution  to  the  kth  time 
fluctuation  in  y,  and  y  represents  the  <|>-t  average.  In  cases  where  the  state  variable  is  a 
function  of  two  independent  spatial  variables,  y=g(<i>,W)>  llte  perturbation  is  given  by 

y(<|>,\|/,t)  =  y(<|>.  V)  +  Xyk(Oak«l>,V) 

k 


Taking  the  first  moment  with  respect  to  <|>  and  t/  yields 

M(,)(t)  =  E^^[y(<|>,\|;,t)]  =  E^[y(<j»,\ir)]4-  y,(t)E^[a,((|>,v)l  +  yj(t)E^[a2(<|>,v)]+.... 

=  JJ  y(<l>.  (<l>.  V)d<l>dV  +  yi(0jj  a,(4>,  v)f4^(<l>.  V)d<|>dV+.... 

♦V  ¥¥ 


=  y+Piyi(0+p2y2(t)+ .  (23) 

The  coefficients  pk  represent  the  spatial  mean  contribution  to  the  kth  tiree  fluctuation  in  y, 
and  y  represents  the  <j)-\jt-t  average.  In  Chapter  6  a  method  for  estimation  of  the  a(z)’s  is 
presented  for  seasonal  orographic  precipitation  and  temperature  fluctuations  in  the  Wasatch 
mountains  of  northern  Utah  (Fan  and  Duffy,  1992). 

In  some  cases  it  will  be  useful  to  relate  the  hypsometric  function  4(<t>)  to 
map  distance  's’  rather  than  elevation.  This  is  done  by  defining  a  characteristic  trajectory  of 
the  catchment,  and  applying  a  change  of  variable.  The  obvious  characteristic  trajectory  is 
the  elevation  of  the  channel,  2(s)  with  projected  distance,  s=(x2+y2)i/2.  To  cover  the  full 
range  of  elevation  in  the  catchment,  the  channel  trajectory  is  extended  to  the  point  of 
maximum  elevation  in  the  catchment  along  the  steepest  path.  Base  level  is  locally  defined 
as  the  outlet  elevation  for  the  catchment,  z=0.  Tlie  specific  area  or  relative  frequency  4(<|)), 
is  now  defined  as  the  fraction  of  total  area  dF^  associated  with  an  increment  of  elevation 
change  d(j>,  which  can  be  related  to  the  projected  distance  ds  along  the  channel  (|>=z(s).  The 
distribution  f<j)(s)  is  given  by 


f4(s)  =  f*[2(s)) 


dz(s) 


ds 


(24) 


The  specific  area  for  the  orthogonal  trajectories  fY(V).  is  also  be  defined  as  in  (24).  In 
terms  of  watershed  geometiy,  f^(s)  or  fy(s)  represent  the  relative  frequency  or  area  per  unit 
change  in  length  along  the  channel.  Thus  f(s)  is  a  measure  of  the  width  or  bounding  shape 
of  the  watershed  relative  to  the  channel.  The  characteristic  shape  of  the  catchment  or 
subcatchment  enclosing  the  longest  channel  is  now  approximated  by  the  pair  of  shape 
functions,  the  characteristic  shape  ^(s)  and  characteristic  elevation  (|>=z(s).  However,  there 
is  no  characteristic  shape  for  hillslopes,  since  a  channel  is  not  enclosed  by  the  hillslope. 
We  use  similarity  arguments  are  used  to  construct  the  characteristic  shape  of  the  hillslope. 

The  perspective  of  landform  'similarity'  was  well  described  in  a  paper  by  Strahler 
(1958),  where  he  stated;  "Systems  of  landforms  involving  the  same  geologic  processes  and 
materials  are  generally  recognized  to  possess  a  considerable  degree  of  similarity.  Indeed 
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the  classification  of  landforms  depends  on  this...”.  Qualitatively,  the  attributes  of  landform 
pattern  and  similarity  imply  an  invariance  with  respect  to  scale  of  measure.  For  example,  in 
small  watersheds  it  is  often  the  case  that  soil  and  colluvium  tend  to  be  thicker  near  the 
stream  than  on  the  ridges,  and  that  the  water  table  tends  to  reflect  a  subdued  image  of  the 
terrain  or  bedrock.  This  pattern  suggests  a  local  form  of  scaling  which  preserves  the  spatial 
relationship  of  field  observations  (soil  thickness)  with  respect  to  natural  boundary 
conditions  (e.g.  ridge  and  stream,  mountain  and  valley).  Duffy,  et  al  (1991)  proposed  a 
local  rescaling  of  topographic  trajectories  to  defme  a  unit  hillslope.  Each  trajectory  is  scaled 
by  its  overall  height  and  distance  from  ridge  to  stream.  The  process  of  rescaling  the  flow 
geometry  and  corresponding  hydrologic  measurements,  leads  to  delineation  of  a  soil 
moisture  and  saturated  storage  envelope,  and  to  an  empirical  analysis  of  the  storage-flux 
dynamics  of  the  hillslope. 

In  general  this  type  of  scaling  is  referred  to  as  self-affine,  as  opposed  to  self¬ 
similar,  since  the  invariance  requires  both  vertical  and  horizontal  scaling  ratios.  For  a  self- 
affine  surface  we  can  say  that  for  any  two  trajectories  z(si)and  z(s2),  the  increments  of 
elevation  and  map  distance  are  proportional  but  with  different  scale  ratios 

dz,  =  X'  dZj  and  ds,  =  X'  ’  dsj 

where  X'  and  X"  represent  scale  ratios  for  the  vertical  and  map  plane  respectively  and 
ds2=dx2+dy2.  This  transformation  is  used  to  preserve  the  relative  position  of  field 
measurements  with  respect  to  boundary  conditions  along  a  path,  and  to  project  sparce  and 
scattered  field  data  onto  a  smooth  curve,  thereby  red.-.ing  the  dimensionality  of  the 
hillslope  from  three  to  two.  Figures  2.3  illustrates  the  concept  of  a  unit  hillslope,  where 
zo(s)  represents  the  rescaled  surface  elevation  trajectories,  zi(s,i)  the  rescaled  location  of 
the  water  table,  and  Z2(s)  the  corresponding  bedrock  position.  Note  that  map  distance  and 
elevation  are  scaled  by  overall  change  for  each  trajectory  As  and  Az  respetively,  s=s'/As 
(()<s<l)  z=z'/Az,  (()<z(s)<l).  The  z(s)'s  are  represented  by  polynomial  shape  functions, 
and  a(t)  and  b(t)  are  length  scales  indicating  the  upslope  penetration  of  the  subsurface 
saturated  region  and  the  surface-saturation  region  respectively.  The  corresponding 
conditional  hypsometric  relation  defines  the  relative  divergence/convergence  of  the 
bounding  trajectories  in  the  map  plane.  In  section  2.4  a  power  law-type  hypsometric 
distribution  is  fit  to  the  Shale  Hills  data. 
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Figure  2.3.  Converging  hillslope  trajectories  defining  a  region  of  local  similarity. 
Similarity  is  tested  by  rescaling  each  trajectory  by  its  overall  height  and  length  to  form  a 
unit  hillslope  zo(s)  trajectory.  The  region  of  similarity  has  trajectories  of  the  same  unit 
shape.  The  corresponding  rescaled  water  table,  bedrock,  partial-area  length  scale,  and  the 
upslope  penetration  of  the  saturated  region  is  indicated. 

2.3  The  Hypsometric  Integral 

A  measure  of  volume  or  fractional  volume  of  land  mass  above  local  base  level  is 

determined  by  integrating  the  hypsometric  distribution  F((j>),  referred  to  as  the  hypsometric 
integral  I, 


A 

I  =  l-jF((}))d<}). 


(25) 


From  the  method  of  derived  distributions  (24),  the  volume  intergation  may  be  applied  to  the 
shape  functions  for  the  land  surface,  water  table  and  bedrock  surface  (zo(s),  zi(s,t),  Z2(s)), 
which  leads  to 


t 

Io=l-jFo(s)ds 

0 

1 

I,(t)  =  l-jF,(s,t 
0 

1 

I2  =  1  -  J  F2(s)ds 


(26) 


(27) 


(28) 


The  differential  volumes  I(rl2  and  I1-I2  represents  the  volume  of  soil  and  the  saturated 
volume  over  the  hydrologic  unit 


2.4  AModelforHypsometry 

In  the  previous  section  probability  arguments  are  related  to  watershed  hypsometry. 
It  was  noted  that  the  problem  of  evaluating  the  local  hypsometric  relation  for  a  hillslope  or 
subcatchment  becomes  an  evaluation  of  the  joint  distribution  of  elevation  and  flow 
trajectories.  To  demonstrate,  a  hypsometric  distribution  is  proposed  of  the  form 
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F(x)=^;  (0^x^l);(a^P^l) 

(a  +  x**) 


(29) 


(30) 


where  x  represents  the  normalized  random  variable,  and  the  distribution  is  defined  for  the 
unit  interval  (O^x^l).  The  specific  area  or  density  f(x)  is  given  by 

<30, 

a  +  x"  (a  +  x”)^ 

The  flexibility  of  (29)  and  (30)  for  describing  watershed  shape  is  illustrated  in  Figure  2.4. 
The  function  captures  a  wide  range  of  hypsometric  shapes  ranging  from  uniformly 
distributed  regions  (rectangular),  to  regions  with  a  large  fraction  of  area  at  high  or  low 
elevation. 


Ill 
Figure  2.4.  Examples  of  the  range  of  shapes  possible  for  the  hypsometric  distribution 
1-F(x)  and  f(x)  from  equation  (29)  and  (30).  From  the  top  row  and  moving  left  to  right  the 
parameters  are:  {{a,p,a})=  {{1.,1.,200.},{1.,1.,.4},{2.,2.,4.},{3.,3.,.01},{7.,7.,.01}, 
{10.,10.,.02}). 


The  topographic  shape  function  z(s)  can  be  represented  as  a  polynomial  in  s 

m 

z(s)  =  XaiS‘ 


(31) 


For  the  Shale  Hills  watershed  2,  z(s)  was  fit  to  a  third  order  polynomial 

z(s)  =  1.027  -  2.738S+ 2.828s^  -  1.128s*  (  32 ) 

The  hypsometric  relations  1-FK<1))  and  1-Fy(\|f)  are  given  in  Figure  2.5.  The  joint 
distribution  F((j>,\|r)  is  shown  in  Figure  2.6. 
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Figure  2.5.  Hypsometric  relationship  for  Shale  Hills  and  the  fitted  model  eq.  (25)  1-F4((|>) 
(a=2,a=2,p=0.4)  and  1-Fy(\|r).  (a=0.009,a=2.4,p=2.3). 


Figure  2.6.  The  joint  distribution  F(<|>,\|f)  for  the  Shale  Hills  watershed  determined  for  the 
hypsometric  relations  shown  in  Figure  2.5. 


The  density  f^((t>)  and  conditional  density  f^(())|A^)  or  width  distribution  for  the  hUlslopes  in 
the  region  adjacent  to  channel  elevation  (0.2^z^.5)  is  shown  in  Figure  2.7. 


2.5  Deriving  the  Dynamical  System  at  the  Hillslope  Scale 

In  this  section  the  hypsometric  distribution  is  used  to  derive  the  dynamical  system 
for  space-averaged  state  variables  at  the  hillslope  scale.  Consider  the  conservation  equation 
for  subsurface  flow  of  moisture  in  the  hillslope  geometry  of  Figure  2.3. 

^  +  Vq’  =  0  (33) 

where  0'  (s,z,t)  =  0(s,z,t)  f^(slA\i/)  represents  the  volumetric  moisture  per  unit  area  in  the 

{s,z}  plane,  0  is  moisture  content  [L^/L^],  f^(slA\)/)  is  the  conditional  density  for  relative 
width  of  the  hillslope,  and  q'(s,z,t)  =  {q’,,q'j}  =  q(s,z,t)  f^(slA\|/)  is  the  volumetric  Darcy 

flow  per  unit  area.  Vertical  integration  of  (33)  is  performed  over  a  column  of  porous  media 


(34) 


The  integration  is  carried  out  in  two  intervals,  below  and  above  the  water  table 
Zj (s)  z  <  zj”  (s,  t)  and  z^ (s,  t)  <  z  ^  Zq  (s),  respectively 

0'  dz  +  ^ dz  =  ^I^f^(slA\|/) -H^:^f^(slAv)  ( 35  ) 

where  11(8,0  and  Y(s,t)  are  vertically  averaged  state  variables  for  the  saturated  and 
unsaturated  depth  of  water  (see  Fig.  2.8).  Applying  Liebnitz  rule  to  the  right  hand  side  of 
(34)  yields 
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*♦  3s  dz 

0q  3q  .  dz7  • .  3z,  ■ ,  3zo  ■ ,  3z^ 

if if 

+q'A.  -qi',;  +q.i„-  -q.'..  <  > 


where 


*1  *0 
q’.,  =  |qldz  and  q'^  =  Jq.dz. 


(37) 


The  final  integration  is  performed  over  the  length  of  the  unit  hillslope  (0<s<l).  Substituting 
for  the  primed  quantities  0'  (s,z,t)  =  0(s,z,t)  f^(slA\jr)  and  q'  (s,z,t)  =  q(s,z,t)  f^(slAv) 

the  unsaturated  and  saturated  storage  terms  become 


1  f 

3t 


3  * 

jil(s,t)  f^(slAv)ds  +  —  jY(s,t)  f^(slA\ir)ds  = 
•«0  **0 


.L  os  .L  3s 


(38a) 

(38b) 


1.0  fO  ««0  ««0 


3z 


3z„ 


3z; 


••0 

1 


(38d) 


+ 1 qil.,  ds -  / Qzl,;  ds  +  J q;L  ds -  J q,l,^  ds 

*4t0  s^O  s>0  (sO 

The  equivalent  depths  of  integrated  soil  moisture  7(s,t)  and  saturated  storage  ii(s,t) 

represent  the  state  variables  of  interest.  Evaluating  the  integrals  in  Ti(s,t)  and  irfs.t)  from 

(38a)  yields  the  conditional  spatial  mean  storage  quantities  for  the  hillslope 

1 

Ti(tlA\|/)  =  J  Ti(s,t)f^(slAv)ds.  ( 39 ) 


tsO 

I 


Y(tlA\|r)=  jY(s,t)f*(slAv)ds  . 


(39) 


tsO 


Equation  (38b)  represents  the  volumetric  boundary  flux  from  the  soil  moisture  and 
saturated  storage  reservoirs  at  the  stream,  q,  +  ql,=q(t).  As  a  first  approximation  we 
assume  that  horizontal  flow  components  in  the  unsaturated  zone  are  small  »  0, 

q' I  »  q' I  and  ~  »  and  thus  the  terms  in  row  (38c)  cancel.  The  terms  in  (38d) 
^  *'  ^  ’  3s  3s 

respectively  are:  spatially  integrated  vertical  flux  evaporation  and  precipitation,  recharge  to 
the  saturated  zone  from  soil  moisture  (r+),  recharge  to  soil  moisture  from  the  saturated  zone 
(r),  and  the  flux  of  water  to  or  from  the  bedrock  (u).  The  dynamical  system  for  the  spatial 
average  storage-flux  over  the  hillslope  is  now  given  by 
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dt  dt 


p-e-r^+r-q-u 


(40) 


Writing  (40)  separately,  combining  surface  inputs,  and  indicating  the  general  dependence 
of  flux  quantities  on  state  variables,  leads  to  a  two-component  state-space  form  of  the 
model 

4^=i(Tl,Y.p,e,t)-r(Ti.y,t)  (41) 

dt 

^=r(Ti,Y,t)-q(Ti,Y.t)  ( 42 ) 

dt 


where  i(n,Y<P>e)  is  surface  infiltration  rate  ,  or  that  flux  which  infiltrates  the  soil  less 
evaporation,  r(T),Y)  is  the  recharge  rate  per  unit  surface  area  to/from  the  water  table,  and 
q(n,Y)  the  volumetric  subsurface  runoff  per  unit  area.  Equations  (41-42)  are  referred  to  as 
the  integral-balance  model  for  the  coupled  unsaturated-saturated  system.  The  state 
variables  are  defmed  locally  in  the  subregion  Ay .  Figure  2.8  illustrates  the  relationship  for 
soil  moisture  6  and  the  vertically  integrated  soil  moisture  depth  y,  and  the  saturated 
thickness  n  for  an  increment  of  hillslope  length  ds  with  surface  elevation  zo(s). 


Figure  2.8  The  soil  moisture  depth  function  Y(s,t)  and  the  saturated  thickness  ii(s,t)  are 
defined  for  a  column  of  porous  material  on  the  hillslope 
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3  SOIL  MOISTURE  NUMERICAL  EXPERIMENTS 


Two  and  three  dimensional  numerical  experiments  were  conducted  based  on  the 
theory  of  unsaturated-saturated  flow  in  a  porous  medium.  The  integrated  storage-flux 
constitutive  relationships  q(il),  r(Y,  t))  and  the  partial-area  length  scale  b(Ti)  were 
determined  as  a  function  of  slope  gradient  (H/L),  slope  curvature/shape,  aspect  ratio  (d/L), 
the  degree  of  topographic  convergence  or  divergence,  soil  type,  and  boundary  conditions. 
Parameterization  is  performed  by  fitting  a  polynomial  to  each  steady-state  relationship  qfq), 
r(Y,  11),  and  sfq).  The  goal  of  the  experiment  is  to  construct  the  functional  form  of  storage- 
flux  and  surface  saturation  relationships  in  the  dynamical  model  as  discussed  in  chapter  2, 
and  to  generally  establish  the  nonlinear  structure  of  spatially  integrated  state  variables. 

3.1  Richards  Equation  and  Solution 

The  partial  differential  equation  describing  this  flow  is  known  as  Richards' 
equation 

C(v)^-V-(K(v)  (Vv+Vz))  (43 ) 

where  y  is  the  soil-water  pressure  head,  C(v)=8e/3y  is  the  specific  moisture  capacity,  0  is 
moisture  content,  K(y)  is  the  nonhysteretic  unsaturated  hydraulic  conductivity,  and  z  is  the 
elevation  head  or  gravity  contribution.  The  Richards  equation  was  solved  by  the  finite 
element  code,  FEMWATER  of  Yeh  (1987),  and  this  reference  gives  details  of  the 
numerical  solution  methodology  based  on  the  Galerkin  finite-element  approximation  as  well 
as  an  input  guide  for  the  program.  An  example  of  element  mesh  design  using  a  few 
quadrilateral  elements  is  illustrated  in  Figure  3.1.  For  most  cases  a  uniform  element  mesh 
distribution  was  used.  The  number  of  elements  for  each  simulation  is  given  in  Table  3.4. 

The  finite-element  formulation  results  in  a  nonlinear  matrix  equation  that  needs  to  be 
solved  by  means  of  iterative  technique.  For  example,  the  following  nonlinear  matrix  is 
obtained  from  the  formulation:  [C(\|r)]  {\|r}  =  {L},  where  [C]  is  the  coefficient  matrix,  \[r  is 
the  dependent  vector  (pressure  head),  and  (L)  is  the  load  vector.  First,  an  initial  guess  of 
for  {t|r)  is  made.  From  this  estimate,  the  nonlinear  matrix  is  linearized  and  solved  by 
Gaussian  elimination.  New  values  for  (v)  are  acquired  using  a  weighted  average  between 
each  solution  and  previous  values  for  This  procedure  is  }  =  {v|;}  +  (1-w)  {\|ri), 
where  j  is  the  iterative  index,  and  w  is  the  iteration  parameter.  Under-relaxation  is 


( 


0<w<l,  exact  relaxation  is  w  =  1,  and  over-relaxation  is  w  >1.  For  all  runs,  the  under- 
relaxation  option  to  Picard's  iteration  scheme  was  the  most  efficient  for  convergence  of  the 
nonlinear  matrix.  The  parameter  w  varied  from  0.1  to  0.5  depending  on  the  case. 


Figure  3.1.  A  quadrilateral  finite-element  mesh. 

In  the  mass-balance  routine,  FEMWATER  calculates  total  moisture  volume  using 
Gaussian  quadrature  integration  technique.  Since  total  moisture  volume  is  a  sum  of 
unsaturated  water  volume  and  saturated  water  volume,  the  integral  unsaturated  water 
storage  was  calculated  as  y=  Ototai  *  tl 

3.2  Flow  Geometry  and  Boundary  Conditions 

The  shapes  and  boundary  conditions  considered  in  the  numerical  experiments  are 
given  in  Figures  3.2  and  3.3.  The  shapes  chosen  here  are  to  some  degree  arbitrary,  but  are 
considered  to  be  representative  of  hillslope  profiles  in  humid-temperate  regions.  A  letter 
(A,  B,  C,  D,  and  E)  is  assigned  to  indicate  each  shape  which  will  be  used  throughout  the 
text.  Table  3.1  defines  the  surface  and  bedrock  shape  functions  that  appear  in  Figure  3.2. 

For  all  hillslope  cases,  the  impermeable  boundary  is  assumed  to  be  B-C.  The  A-B 
and  C-D  are  symmetry  boundary  conditions  such  that  no  flow  crosses  these  boundaries. 
The  constant  pressure  head  condition  at  the  point  D  represents  a  perennial  stream.  The 
infiltration  and  seepage  boundary  conditions  coexist  on  the  soil  surface  (A-D),  and  the 
position  of  the  seepage  face  is  not  known  a  priori.  The  boundary  condition  is  tested  during 
each  solution  (see  Yeh,  1987). 
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For  the  horizontal  stream-aquifer  system  (shape  E),  the  no-flow  boundary  condition 
was  imposed  on  the  boundaries  through  A-B-C  and  D-E.  A  fully  penetrating  stream  was 
assumed  on  the  boundary  C-D,  and  a  constant  pressure  head  represents  the  stream.  On  the  • 

land  surface  A-E,  a  variable  boundary  condition  was  specified  to  handle  the  infiltration  rate. 

A  variable  condition  means  a  mixed  boundary  condition  for  Neuman  (flux)  type  and 
Dirichlet  type.  For  each  iteration,  either  Neumann  type  or  Dirichlet  type  is  determined  by 
flux  balance  consideration  between  rainfall  flux  and  Darcy  flux. 


I  i  m  I  i 


On  A-B-C-D;  q  n=0 
OnA-D;  q  n=P  n  or  \|r(x,t)=0 
At  the  point  D;  vr  =  0 


Figure  3.2.  Hillslope  shapes  and  boundary  conditions  (p  =  rainfall  flux  vector,  q  = 
Darcy  flux  vector,  and  n  =  unit  normal  vector). 
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Shane  E 


I  I  I  I 


A  E 


X 


Figure  3.3.  A  horizontal  aquifer  shape  and  boundary  conditions  ( p  =  rainfall  flux 
vector,  q  =  Darcy  flux  vector,  and  n  =  unit  normal  vector). 

Table  3.1,  Surface  and  bedrock  shape  functions. 
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Table  3. 1  (cont.) 


3.3  Soil  Texture  and  Hydraulic  Properties 


In  general,  soil  hydraulic  properties  vary  with  moisture,  soil  texture  and  soil 
structure.  The  numerical  experiments  were  conducted  for  four  soils  (Guelph  Loam, 
Plainfield  Sand,  Yolo  Clay,  and  Touchet  Silt  Loam  G.E.3).  Soil  hydraulic  properties  for 
these  soils  reported  in  Elrick  et  al.  (1990)  and  van  Genuchten  (1980)  are  reproduced 
below.  The  unsaturated  hydraulic  conductivity  K(v)  for  Guelph  Loam,  Plainfield  Sand, 
and  Yolo  Clay,  is  given  by 


K(v)=K,exp(a\|/) 
and  for  Touchet  Silt  Loam  G.E.3 

ll-(a\|i)"-'ll+(ai|()"r")^ 

— [i+(av)"r'' — 

where  Kg  is  the  saturated  hydraulic  conductivity,  and  a  ,  a,  n,  m=l-l/n  are  fitting 
parameters.  Hie  relationship  between  soil  moisture  content  and  matric  potential  is  given  for 
Guelph  Loam,  Plainfield  Sand,  and  Yolo  Clay,  as 

V<0 


e=  ,P)P; _ +D 

and  for  Touchet  Silt  Loam  G.E.3,  as 


e=e-+(0g-0,)0,  0=  - - — - 


\|;<0 


Soil  water  capacity  relations  C(v)  are  obtained  by  differentiating  the  water  retention  curve. 
The  for  Guelph  Loam,  Plainfield  Sand,  and  Yolo  Clay  have  the  form 
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(Pj-» 


C(w)=ii=P!MaM _ 

and  the  Touchet  Silt  Loam  G.E.3  is 

C(\|;)=— -^1.)  Qt/m  jm 

dvjf  1-m 


V<0 


V  <0 


Table  3.2  summarizes  parameter  values  for  each  relationship.  Figure  3.4  illustrates  soil 
moisture  retention  and  relative  hydraulic  conductivity  curves  as  a  function  of  pressure  head. 


Table  3.2.  The  oarameter  values 


soil  type 

Kg  ( m/hr ) 

6s 

a  ( 1/m) 

P] 

P2 

P? 

P4 

Guelph  Loam 

0.013212 

0.523 

3.36 

0.243 

0.421 

2.0 

0.28 

Plainfield  Sand 

0.12384 

0.477 

13.06 

0.377 

0.00154 

4.0 

0.1 

Yolo  Clay 

0.0004428 

0.499 

1.88 

0.259 

0.609 

1.5 

0.24 

soil  type 

Kg  ( m/hr ) 

6s 

6r 

a 

n 

Touchet  Silt  Loam  G.E.3 

0.12625 

0.469 

0.19 

0.243 

0.421 

Figure  3.4.  Water  retention  function  0(\|/)  and  relative  conductivity  function  K(tjr)/Ks  for 
Guelph  Loam  (1),  Plainfield  Sand  (2),  Yolo  Clay  (3),  and  Touchet  Silt  Loam  G.E.3  (4). 

•  3.4  Equilibrium  Experiments 

Equilibrium  numerical  experiments  were  conducted  to  evaluate  the  steady-state 
(constant  forcing)  constitutive  relations  for  the  spatially-integrated  state  variables  (q  and  y) 
and  the  relation  to  system  fluxes.  The  partial  area  length  scale,  or  region  of  saturation 
®  overland  flow  was  also  determined.  The  experiments  were  comprehensive  in  that  for  each 
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geometry  and  parameterization,  constant  rainfall  rates  were  specified  to  cover  all  possible 
cases  of  soil  moisture  and  saturated  storage,  from  very  dry  and  entirely  unsaturated  to 
complete  saturation.  Table  3..3  lists  physical  information  and  Table  3.4  describes 
computational  information  for  20  simulated  cases. 


Note:  L  =  slope  length  =  100  m  and  H  =  slope  height  NA  =  not  applicable. 
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Table  3.4.  Computational  information  for  steady-state  cases. 


case 

nodes 

■wei 

dx/dz 

w 

1 

2211 

0.5 

2 

0.25 

2 

M 

M 

II 

H 

3 

M 

N 

H 

M 

wm 

H 

H 

n 

0.15 

5 

II 

M 

II 

0.25 

6 

H 

II 

n 

N 

7 

N 

n 

It 

n 

8 

tl 

H 

N 

H 

9 

3276 

M 

N 

0.1 

10 

3861 

H 

II 

II 

11 

4221 

0.5 

4 

0.15 

12 

2211 

M 

2 

0.25 

13 

4221 

H 

4 

0.15 

14 

4141 

1.0 

2 

0.1 

15 

II 

H 

II 

0.05 

16 

2211 

0.5 

2 

0.25 

17 

ri 

II 

II 

II 

18 

II 

H 

II 

II 

19 

II 

•I 

II 

II 

20 

II 

II 

.  If 

II 

Note:  dx  =  horizontal  length  of  an  element,  dz  =  vertical  length  of  an  element,  w  is  the 
iteration  parameter  for  solving  a  nonlinear  matrix  equation. 

The  cross-section  shape  and  symbols  in  Table  3.3  refers  to  Figures  3.2  and  3.3. 

Changes  in  slope  gradient  alone  (H/L)  were  considered  in  cases  1-4,  while  other 
parameters  were  fixed.  For  cases  5  and  6,  the  aspect  ratio  (d/L)  was  varied  while  the 
remaining  parameters  were  the  same  as  case  1.  Cases  7  and  8  examine  the  effect  of 
different  hilllslope  shapes  (linear  and  concave).  Cases  9  and  10  investigate  the  role  of 
storage  capacity  on  the  flux-storage  relationship  by  changing  the  depth  to  bedrock.  The 
effect  of  soil  type  was  examined  in  cases  1 1-13,  and  cases  14  and  15  simulate  a  horizontal 
hillslope  system  with  different  soil  types. 


Cases  16*19  show  the  effect  of  convergence  or  divergence  in  plan  or  map  view, 
while  the  remaining  parameters  follow  case  1.  Figure  3.5  shows  the  converging  and 
diverging  sections  for  these  cases.  The  degree  of  convergence  or  divergence  was 
controlled  by  the  parameter  p/L,  that  is,  the  ratio  of  the  radius  of  curvature  in  the  plan  view 
(p)  to  the  hillslope  length  (L).  Since  a  circular  converging  or  diverging  domain  is 
considered,  p  is  simply  the  radius  of  a  circle.  For  other  planametric  shapes,  p  can  be 
defined  as  the  inverse  of  the  magnitude  of  curvature  vector  and  depends  on  position. 
Notice  that  p/L  approaches  infinity  for  a  parallel  hillslope,  or  the  hillslope  length  is  small 
relative  to  p.  The  ratio  p/L=l  is  the  other  limiting  case,  and  the  range  is  1<  p/L  <«.  For 
each  plan  shape,  two  cases  were  examined,  p/L=l  and  p/L=2.  Solutions  were  obtained 
using  axisymmetric  coordinates  in  FEMWATER. 

Case  20  tests  the  effect  of  the  stream  boundary  condition.  Instead  of  a  constant  head 
condition,  a  seepage  face  was  defined  as  shown  in  Figure  3.6.  All  other  parameters 
remained  the  same  as  case  1.  Case  1  represents  the  base  case  for  comparison. 


Figure  3.5.  Plan  view  of  hollows  and  spurs  (F-F  =  stream  and  G-G'  =  divide). 


On  A-B-C;  q .  n=0 

OnA-D-C;  q  n=P  n  or  \|/(x,t)=0 
Figure  3.6.  Boundary  conditions  for  case  20. 

3.4.1  Saturated  Storage-Runoff  Relation 

Figure  3.7  summarizes  the  steady-state  subsurface  outflow  or  baseflow  versus 
saturated  storage,  q(T|)  for  a  complete  range  of  (constant)  precipitation  rate  p,  and  all  cases 
in  Table  3.4.  Each  point  represents  a  single  steady-state  solution  to  the  Richards  equation 
for  T|  and  q  with  specified  p.  A  third-degree  polynomial  was  fit  to  the  calculated  values  of 
q(Ti).  Polynomial  coefficients  for  all  cases  are  summarized  in  Appendix  A. 

Note  that  for  steep  slopes,  the  shape  of  q(Ti)  tends  to  become  more  convex  (case  4). 
The  magnitude  of  q(Ti)  is  similar  for  cases  1,  2,  and  3.  For  dry  to  moderately  wet 
conditions,  a  linear  q(Ti)  can  be  defined  for  convex-concave  hillslopes,  at  least  for  the  range 
of  conditions  examined  here.  For  steep  slopes  however,  the  magnitude  of  q  is  a  nonlinear 
function  of  q. 

When  the  aspect  ratio  of  the  hillslope  (d/L=soil  depth/slope  length)  is  increased, 
almost  no  change  in  the  shape  of  q(q)  is  observed.  However,  the  magnitude  of  q(q) 
increases  rapidly  with  aspect  ratio  (cases  1, 5,  and  6). 
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case  19 
p/L=2 


0.2  0.7  1.2 


X 


q  (mm/hr) 


case  20 


0.7  1.3 


Ti(m) 


Figure  3.7.  Steady-state  q(Ti)  relationships.  Dots  were  determined  from  steady-state 
solutions  to  Richards  equation  for  a  constant  rainfall  rate.  The  solid  curve  is  the  third 
degree  fitted  polynomial.  Soil  type  assumed  is  soil  2  =  Plainfield  Sand,  soil  3  =  Yolo 
Clay,  and  soil  4  =  Touchet  Silt  Loam  G.E.3.  Guelph  Loam  soil  is  assumed  unless 

specified.  The  fitted  polynomial  is  the  third-order:  q(ri)  =  ao  +  ai  t]  +  a2  ri^  +  33  T|^. 
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Hillslope  curvature  (cases  1,  7,  8)  substantially  alters  the  shape  of  q.  The  linear 
slope  (case  7)  exhibits  a  concave  q-T]  relationship,  while  cases  1  and  7  have  slightly  convex 
to  concave  q-n.  For  the  linear  slope,  q  approaches  saturation  at  a  lower  value  of 
(q  >0.7m),  as  compared  to  the  convex-concave  or  concave  slopes,  indicating  the  greater 
efficiency  of  linear  slopes  in  generating  runoff. 

When  depth  to  bedrock  increases  and  storage  capacity  increases,  a  concave  relation 
for  q(T])  emerges  (cases  9  and  10).  The  magnitude  of  q  also  increases.  This  relationship 
could  be  quite  useful  in  field  studies  where  the  storage  geometry  is  unknown.  Ando  et  al. 
(1983)  found  a  similar  relation  q(q)  in  a  forested  hillslope  watershed,  and  fitted  a  quadratic 
polynomial  to  saturated  storage  and  discharge.  As  the  aspect  ratio  d/L  increases  (case  10), 
the  shape  of  the  relationship  displays  somewhat  less  concavity  and  the  range  of  q  increases 
slightly. 

The  effect  of  soil  type  on  q(T})  was  examined  for  hillslopes  in  cases  11-13  and  for 
horizontal  flow  in  cases  14-15.  For  all  cases  soil  type  did  not  affect  the  basic  form  of  q(Ti) 
(cases  1, 11, 12,  13,  14  and  15).  However,  variations  in  soil  type  have  a  large  impact  on 
the  magnitude  of  q.  The  increase  in  the  magnitude  of  q  was  found  to  be  nearly  proportional 
to  the  increase  in  saturated  hydraulic  conductivity  for  a  given  aquifer  shape. 

Cases  16-19  illustrate  the  effect  of  topographic  convergence  and  divergence  on 
q(Ti).  As  noted  previously,  topographic  convergence  or  divergence  in  plan  (Figure  3.5) 
was  simulated  for  an  axisymmetric  coordinate  system.  From  radial  symmetry,  the  soil 
moisture  gradient  is  constant  along  radii  from  the  center  of  convergence.  In  order  to  fairly 
compare  the  response  with  a  parallel  case,  the  total  hillslope  volume  for  the  hollow  or  spur 
geometry  was  made  the  same  as  the  parallel  slope.  As  might  be  expected  for  radial 
symmetry,  the  simulations  demonstrate  the  similarity  in  the  shape  of  q(Ti),  independent  of 
the  overall  angle  of  topographic  convergence  or  divergence.  For  the  same  1),  as  p/L  ratio 
increases,  q  decreases  for  hollow  cases,  and  increases  for  spur  slopes.  Also,  as  the  ratio 
p/L  increases,  the  relationship  q(T|)  for  hollow  or  spur  approaches  the  limiting  case  of  a 
parallel  slope  (p/L  = «»)  as  expected. 

Case  20  displays  the  effect  of  boundary  condition  at  the  stream.  Compared  to  case 
1,  the  shape  of  q(T|)  in  case  20  is  quite  different  for  dry  and  very  wet  moisture  conditions, 
but  for  the  medium  wetness  range,  the  two  cases  exhibit  the  same  linear  relationship  q(T|). 
The  amplitude  of  subsurface  outflow  in  case  20  is  slightly  larger  than  case  1. 


With  the  same  series  of  steady-state  simulations,  recharge  rates  were  estimated  and 
parameterized  in  terms  of  y  and  t].  Figure  3.8  plots  these  results  for  r(Y,  t^)  and  its 
projection  onto  three  orthogonal  planes:  r-y  plane,  r-T)  plane,  and  y-i]  plane.  The  a  cubic 
polynomial  was  fit  to  r(y,  t))  and  the  results  are  given  in  the  Appendix. 

Figure  3.8  shows  that  r  decreases  with  an  increase  in  y.  The  soil  moisture  y  generally 
decreases  with  an  increase  in  t).  Both  y-r  and  y-T)  exhibit  an  inverse  relationship  over  most 
of  the  range  of  rainfall.  Recharge  r  is  directly  proportional  to  saturated  storage  t].  This 
can  be  interpreted  to  mean  that  for  an  incremental  increase  n  precipitation,  the  saturated 
zone  will  increase  at  the  expense  of  unsaturated  moisture  storage.  The  implication  is  that  y- 
r\  are  competitive  for  the  overall  storage  volume.  This  is  an  important  structure  with  major 


implications  to  modeling  a  two-state  system.  For  very  low  moisture  conditions  however,, 
the  relationships  y-r  and  y-il  are  direct  rather  than  inverse.  Under  very  dry  conditions,  we 
expect  both  soil  moisture  and  saturated  storage  to  increase.  At  some  critical  value  of  y-T)  a 
threshold  is  reached,  and  the  relationship  becomes  inverse  as  p  increases.  This  effect  may 
be  extremely  important  for  explaining  the  transition  from  drought  to  wetter  conditions  or 
vice-  versa,  as  the  constitutive  relation  apparently  changes  form. 

The  recharge  relationship  r(y,T})  exhibits  only  modest  changes  in  shape  and 
magnitude  when  the  geometry  and  physical  attributes  of  hillslopes  are  varied.  Chapter  5 
will  take  up  this  discussion  again  in  the  context  of  model  development  and 
parameterization. 
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Figure  3.8.  Steady-state  r(7,T|)  relationships  for  steady-state  solutions  to  the  Richards 
equation  (Lee,  1993).  The  dots  represent  the  steady- state  solutions  for  a  given  precipitation 


Figure  3.8  (cont.) 

3.4.3  Paniai  Area-Saturated  Storage  Relation 

As  defined  earlier,  b  is  the  length  or  area  fraction  of  surface  saturation  for  the 
hillslope.  The  length  scale  b  was  parameterized  as  a  function  of  the  saturated  storage  (t|). 
In  the  numerical  simulations,  b  was  estimated  as  the  distance  over  which  atmospheric 
pressure  is  maintained,  as  illustrated  in  Figure  3.9. 

Figure  3.10  shows  the  relationship  b-ii  obtained  from  the  numerical  simulations. 
The  parameterization  was  performed  by  fitting  a  third-degree  polynomial  in  Tj,  and  results 
are  given  in  Appendix  C. 

The  effect  of  relief  H/L  on  b(Ti)  is  apparently  small  (cases  1, 3  and  5),  with  nearly 
the  same  shape  and  magnitude  for  all  cases.  The  shape  of  b(T|)  between  different  aspect 
ratios  is  similar  despite  the  difference  in  the  scale  of  ii  (case  1,  case  5,  and  case  6).  The 
slope  curvature  has  a  relatively  large  effect  on  the  shape  of  b(ii)  (case  1,  case  7,  and  case 
8).  In  case  7  (straight  slope),  b(o)  is  concave,  with  b  increasing  rapidly  with  t|.  Over 
most  of  the  range  of  t),  the  magnitude  of  b  for  a  linear  slope  is  smaller  than  for  concave  or 
convex-concave  slopes. 

Comparing  the  effect  of  soil  thickness  of  cases  9  and  10,  we  see  that  the  shape  of 
b(Ti)  is  similar.  However,  the  penetration  distance  (b)  is  much  larger  than  is  the  case  for 
shallow  soils.  Soil  type  has  very  little  impact  on  the  shape  and  magnitude  of  b(Ti),  as 
observed  in  cases  11, 12  and  13. 
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The  effect  of  topographic  convergence  and  divergence  on  the  shape  and  scale  of  b(Ti)  is 
shown  in  cases  16-19.  For  converging  slopes  (hollow)  with  p/L  =1  and  2  (cases  16  and 
17),  the  shape  and  amplitude  of  b(Ti)  is  quite  close  to  the  parallel  slope  (case  1).  For 
diverging  slopes  (spur)  the  magnitude  of  q  is  reduced  near  the  channel. 

Case  20  (boundary  condition  effect)  is  similar  to  case  7  (linear  slope),  where  growth 
of  the  zone  of  saturation  is  limited  until  saturated  storage  is  quite  high,  ti>0.4.  The  shape 
of  b(q)  is  concave  for  both  cases.  By  comparison  with  case  1,  where  the  boundary 
condition  is  constant  head  or  stream  stage,  the  partial  area  b  will  increase  even  at  low  values 
of  saturated  storage. 
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Figure  3.10.  Steady-state  b(ii)  relationships.  Dots  were  detennined  from  steady-state 
solutions  to  Richards  equation.  The  solid  curve  is  the  third  degree  fitted  polynomial.  Soil 
type  assumed  is  soil  2  =  Plain  Held  Sand,  soil  3  =  Yolk  Clay,  and  soil  4  =  Touchet  Silt 
Loam  G.E.3.  Guelph  Loam  soil  is  assumed  unless  specified.  The  fitted  polynomial  is  the 

third-order:  bfq)  =  ao  +  ai  q  +  a2  +  a3  r\^. 
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^■4.4  Summary  of  Eauilibrium  Experiments 

The  steady-state  storage-flux  relations  q(tl).  r(Y,  t|),  and  b(ti)  were  established 
through  numerical  experiments  with  Richards  equation  for  wide  range  of  hillslope 
geometry  and  hydraulic  properties.  The  shape  of  the  steady-state  q(n)  relation  is  mainly 
affected  by  the  degree  of  slope  gradient,  hillslope  shape  (curvature),  and  aquifer  storage 
capacity.  On  the  other  hand,  the  effects  of  aspect  ratio,  soil  type,  and  topographic 
convergence  or  divergence  on  the  shape  of  q(q)  are  small  for  a  given  hillslope  geometry. 
The  shapes  for  q(q)  show  a  combination  of  linear,  convex,  and  concave  components. 
Moreover,  the  magnitude  of  qCq)  is  strongly  influenced  by  the  saturated  hydraulic 
conductivity,  and  the  capacity  of  saturated  water  storage  or  aspect  ratio. 

Steady-state  recharge  rate  was  examined  in  terms  of  both  state  variables  y  and  q.  The 
relationship  r(Y,  q)  is  determined  by  the  geometric  attributes  of  slope,  hillslope  shape,  soil 
type,  storage  capacity,  and  boundary  condition  near  the  stream.  The  relationships  y-r  and 
y-q  are  inverse  for  wet  soil  moisture  conditions.  Apparently  y-q  are  competitive  for  the 
available  storage  volume.  However,  y-r  and  y-q  have  a  direct  relationship  under  dry 
conditions.  The  implications  of  this  threshold  and  transition  are  critical  to  modeling  and  are 
discussed  further  in  chapter  5. 

In  general,  the  surface  saturation  area  increases  with  the  increase  in  the  saturated 
water  storage.  For  the  same  slope  shape,  the  effect  of  aspect  ratio,  soil  type,  and  storage 
capacity  on  the  shape  of  steady-state  b(q)  is  small,  whereas  the  shape  of  b(q)  changes  with 
the  degree  of  slope  gradient  and  topographic  convergence  or  divergence.  The  shape  of 
b(q)  is  also  affected  by  hillslope  shape  and  the  boundary  condition  at  the  stream.  The 
range  of  b  is  always  bounded  between  0  and  1,  while  the  range  of  q  depends  on  porosity 
and  aquifer  storage  capacity. 

The  motivation  for  this  reasonably  comprehensive  series  of  steady-state  numerical 
experiments  is  to  improve  the  conceptual  basis  for  interpreting  storage-flux  relationships 
observed  in  the  field.  The  nature  of  dynamic  forcing  is  examined  in  the  next  chapter. 

3.5  Transient  Experiments 

The  steady-state  storage-flux  relationships  examined  in  section  3.4  provide  insight 
into  the  equilibrium  patterns  of  hydrologic  variables  as  a  function  of  geometry  and 
hydraulic  parameters.  (Dne  expects  however,  that  the  general  storage-flux  relationships  will 
be  time  dependent.  Numerical  experiments  in  this  chapter  investigate  the  dynamics  of 
subsurface  storage-flux  subject  to  variations  in  the  input  parameters.  Special  attention  is 
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given  to  determining  the  dynamic  relationships  for  storage  y-Tl.  runoff  q(T|,t)  and  recharge 
r(Y.Tl,t).  The  dynamic  experiments  are  by  no  means  comprehensive,  and  additional 
experiments  are  planned. 

The  dynamic  experiments  are  conducted  for  two  geometries,  horizontal  flow  and 
convex-concave  hillslopes.  The  state  variables  are  presented  in  two  ways.  The  first  is 
simply  the  time  series  of  the  response  and  the  second  is  phase  plane  or  phase  portrait  of  the 
response.  The  phase  portrait  represents  a  graph  of  the  trajectories  of  the  state  variables  y-t| 
with  time  implicit.  Chapter  S  provides  a  brief  discussion  of  the  phase  portrait  as  basic  tool 
of  dynamical  systems  analysis.  The  constitutive  relations  q(Ti)  and  r(Y,Ti)  are  also  plotted 
against  the  state  variable  >vith  time  implicit.  The  structure  of  y-T)  and  q-T|  trajectories  is  first 
examined  using  a  simple  step  input  from  initially  steady-state  conditions.  The  effects  of 
finite-duration  rainfall  intensity,  hillslope  shape,  soil  type,  and  boundary  conditions  are 
investigated  to  understand  their  role  on  storage-flux  dynamics.  Finally,  the  Y*'n  *tnd  q-T) 
trajectories  are  studied  using  a  periodic  forcing  (sinusoid),  and  the  response  characteristics 
of  storage  and  fluxes  are  analyzed  and  interpretations  are  made.  For  the  dynamic 

experiments  the  recharge  rate  or  flux  to  or  from  y  and  T|,  was  estimated  from  the  numerical 

dn 

solution  using  the  following  form  of  mass  balance  for  the  saturated  zone,  r=-^+q. 
Given  discrete  values  for  i),  and  qi,  r,  is  calculated  explicitly  r,^j  =  -  ti,  )  /  dt  +  q,^, . 

3.5.1  Step  Input  Response 

For  the  step  input  shown.  Table  3.5  describes  the  particular  information  about 
simulation  shape  and  input  characteristics  for  each  of  the  four  cases  examined  in  this 
section,  and  Figure  3.11  illustrates  the  input  form  and  initial  conditions.  The  simulation 
shape  and  boundary  conditions  for  the  horizontal  flow  case  are  shown  in  Figure  3.3.  For 
the  hillslope  cases,  shape- A-and  the  boundary  conditions  of  Figure  3.2  were  chosen.  The 
geometric  and  physical  parameter  values  used  in  case  1  and  case  2  in  Table  3.5  are  the  same 
as  case  14  in  Table  3.3,  and  the  parameter  values  for  cases  3  and  4  in  Table  3.5  are  the 
same  as  case  1  in  Table  3.3.  The  initial  condition  for  each  experiment  was  based  on  the 
steady-state  solution  obtained  for  rainfall  intensity  fo. 

Simulation  results  are  presented  in  Figure  3.12,  and  the  Y-fl  phase  trajectories  are 
found  to  closely  follow  the  steady-state  results,  which  are  also  shown.  Note  that  for  all 
cases,  the  yfl  phase  trajectories  have  an  inverse  relation  except  at  very  low-storage 
conditions. 
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Figure  3.1 1 .  The  rising  step  (+)  and  falling  step  (-)  for  rainfall. 
Table  3.S.  Aquifer  shape  and  forcing  feature  of  step  input  cases. 


case 

_ type _ 

foflCs 

a/Ks 

1 

0.001 

0.03 

horizontal  rectangular 

2 

0.03 

0.001 

II 

3 

0.001 

0.02 

hillslope 

MM 

0.02 

0.001 

II 

note:  Ks  =  0.013212  m/hr  =  saturated  hydraulic  conductivity  of  Guelph  Loam  soil. 

Figure  3.13  illustrates  the  transient  q-Ti  relationship.  The  transient  relationship  for  q- 
y\  also  follows  the  steady-state  result  fairly  well. 

Figure  3.14  represents  the  time  series  of  infiltration,  recharge,  and  subsurface 
outflow  for  the  horizontal  flow  system.  Case  1  represents  the  (+)  step  change  in  input,  and 
case  2  represents  the  (-)  or  recession.  It  was  later  determined  that  the  rapid  fluctuation  in 
recharge  were  the  due  to  the  explicit  equation  used.  The  fluctuations  dissipate,  but  do 
somewhat  confuse  the  interpretation  of  the  recharge  response.  In  case  1  with  a  rising  step, 
there  is  no  change  in  the  integral  recharge  rate  until  about  1440  hours,  after  which  the 
recharge  rate  increases  rapidly  reaching  a  peak  nearly  2400  hours  after  the  start  of  the 
event. 

During  the  recession,  recharge  decays  slowly  and  converges  to  the  equilibrium 
infiltration  rate.  It  is  interesting  to  note  that  recharge  rate  exceeds  infiltration  during  most  of 
rising  phase.  That  is,  the  net  recharge  rate  has  been  enhanced  by  existing  soil  moisture  in 
unsaturated  zone.  This  demonstrates  the  importance  of  soil  moisture  in  controlling  the 
recharge  dynamics.  Similar  to  recharge,  subsurface  outflow  exhibits  a  delay  at  early  time, 
and  then  increases  to  steady  state. 
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Figure  3.14.  Infiltration  (i),  recharge  (r),  and  subsurface  outflow  (q)  responses  to  cases  1 
and  2.  Note:  the  oscillations  in  r  are  an  artifact  of  the  method  used  for  estimation  . 


For  case  2,  with  a  falling  step,  recharge  quickly  decreases  during  the  event,  and  the 
sign  of  recharge  changes  from  positive  to  negative.  Negative  recharge  (loss  to  saturated 
zone)  is  maintained  during  the  recession  to  steady  conditions.  Negative  recharge  represents 
a  loss  from  saturated  to  unsaturated  storage,  where  the  water  table  is  relaxing  while  leaving 
behind  soil  moisture.  The  subsurface  outflow  response  in  case  2  shows  a  smooth 
relaxation,  while  the  recharge  response  exhibits  a  slow  response  with  numerical  oscillations 
as  in  case  1. 

Figure  3.15  shows  the  time  series  for  infiltration,  recharge,  and  subsurface 
outflow,  and  the  surface  saturation  area  for  the  hillslope  response  resulting  from  a  rising 
and  falling  step  input.  In  case  3  (rising  step)  the  infiltration  rate  increases  with  time  until  it 
reaches  infiltration  capacity,  and  then  decreases  smoothly  with  time.  The  reduction  in 
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inHltration  rate  is  caused  by  the  expansion  of  the  surface  saturation  area  where  rainfall  is 
rejected.  The  recharge  rate  displays  a  bimodal  distribution  in  time,  with  a  maximum  value 
at  around  1000  hours  for  this  soil.  The  q  and  b  response  both  increase  with  time.  • 

In  case  4  (hillslope-falling  step),  recharge  changes  sign  from  positive  to  negative 
similar  to  the  horizontal  flow  case,  and  remains  negative  during  recession.  The  q  and  b 
response  is  a  simple  decay.  The  investigation  next  examines  the  dynamics  of  storage  and 
flux  subject  to  time  variations  in  rainfall  pattern,  and  the  role  of  slope  shape,  soil  type,  and 
boundary  conditions.  ^ 


1000  2000 
time  (hours) 


Figure  3.15.  Infiltration  (i),  recharge  (r),  subsurface  outflow  (q),  and  partial  area 
response  (b)  to  cases  3  and  4. 
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3.5.2  Finite  Duration  Step  Innut 

For  a  discrete  rainfall  duration  as  shown  in  Figure  3.16,  three  different  cases  were 
investigated  by  changing  the  intensity  or  duration  of  rainfall  Table  3.6  describes  input  data 
for  the  intensity  (a)  and  duration  (S)  for  each  case.  The  hillslope  shape-A  and  boundary 
condition  of  Figure  3.2  were  applied  to  all  cases,  the  soil  was  a  Guelph  Loam,  H/L=0.2S, 
d/Ls0.025,  and  LslOO  m  were  assigned.  The  same  parameters  as  case  1  in  Table  3.3  were 
applied  to  all  pulse  input  cases.  To  obtain  a  starting  initial  condition,  a  steady-state 
simulation  was  made  with  rainfall  intensity  fo.  After  the  storm  of  duration  5  and  uniform 
rainfall  intensity  of  a,  the  experiment  was  continued  until  the  solution  returned  to  the 
original  equilibrium  (~250  days). 


Figure  3.16.  Finite  duration  rainfall  event. 
Table  3.6.  Pulse  input  cases. 


case 

fo/Ks 

o/Ks 

6  (hours) 

5 

0.002 

0.2 

48 

6 

M 

0.4 

M 

n 

If 

0.2 

24 

note:  Ks  =  0.013212  m/hr  =  saturated  hydraulic  conductivity  of  Guelph  Loam  soil. 


Figure  3.17  shows  the  transient  response  of  q-T|  for  all  3  cases.  Clearly,  the 
intensity  and  duration  of  the  rainfall  event  effect  the  response.  The  departure  of  the 
transient  relationship  from  steady-state  is  also  indicated.  This  departure  is  large  when  either 
the  intensity  or  the  duration  increases.  The  q-T|  trajectory  exhibits  hysteresis. 

Figure  3.18  illustrates  the  phase  trajectories  for  y-Tj,  which  exhibit  clockwise 
looping.  The  looping  size  of  y-il  trajectories  varies  with  rainfall  intensity  or  duration.  The 
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looping  pattern  itself  is  similar  in  shape  for  all  cases.  Figure  3.19  shows  time  series  for  i, 
r,  q,  and  b  for  the  pulse  input  simulations.  The  infiltration  time  series  for  all  3  cases  is 
similar,  increasing  smoothly  until  infiltration  capacity  is  reached  at  7  to  8  hours,  and  then 
decreasing  with  the  expansion  of  surface  saturation  area. 

Apparently  subsurface  runoff  q  is  characterized  by  two  peaks.  The  initial  peak  for 
all  3  cases  occurs  at  the  end  of  the  input  event  The  first  peak  magnitude  varies  with 
rainfall  volume.  For  case  6,  with  the  largest  event  volume,  the  magnitude  of  the  delayed 
peak  is  greater  than  the  first  peak.  For  case  7  with  shorter  duration/volume,  the  delayed 
peak  is  much  smaller.  Clearly  the  magnitude  and  timing  of  the  delayed  peak  (more  than  1 
month  delay  after  the  storm)  depends  on  the  rainfall  characteristics,  but  it  is  also  a  function 
of  the  low  saturated  hydraulic  conductivity  of  Guelph  Loam  soil  used  in  the  simulations. 
Additional  simulations  for  the  remaining  soil  types  were  conducted.  Remarkably  the  timing 
of  the  second  peak  remained  >10(X)  hours,  but  the  magnitude  of  the  second  peak  decreased 
for  more  permeable  soils.  This  large  delay  in  the  second  peak  for  these  simulations 
illustrates  an  important  nonlinear  response  of  Richards  equation. 

Changes  in  the  recharge  response  of  Figure  3.19  are  most  dramatic  during  the  rise 
and  fall  of  the  pulse  input.  Again,  small  oscillations  are  due  to  the  estimation  method.  The 
recharge  rate  in  case  6  exhibits  a  distinct  pattern  compared  to  other  cases  in  that  the  peak 
recharge  occurs  slightly  after  the  storm  event.  Note  that  the  peak  recharge  rate  exceeds 
infiltration  capacity,  as  was  noted  for  the  step  input  cases. 

The  dynamics  of  the  surface  saturation  area  b(t)  exhibit  nearly  the  same  response  as 
subsurface  outflow.  The  timing  of  the  initial  and  delayed  peak  is  identical  for  q(t)  and  b(t). 
Again,  the  expansion  or  contraction  of  b(t)  largely  depends  on  rainfall  intensity  and 
duration  if  other  conditions  on  the  hillslope  are  fixed. 
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Figure  3.17.  Trajectories  for  q  and  q  for  a  pulse  input.  Numbers  denote  the  position  of 
path  in  days  from  the  start  of  simulation.  Dots  refer  to  steady-state  simulations  from  case  1 
in  Table  3.3. 
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Figure  3.18.  Phase  trajectories  for  y-T).  Numbers  denote  the  position  of  path  in  days  from 
the  start  of  simulation.  Dots  refer  to  steady-state  simulations  from  case  1  in  Table  3.3. 
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Figure  3.19.  Infiltration  (i),  recharge  (r),  subsurface  outflow  (q),  and  fraction  of  surface 
saturation  area  (s)  responses  for  cases  S-7. 
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^  In  this  section  the  numerical  experiments  examine  periodic  variations  in  rainfall. 

Moisture  conditions  fluctuate  from  relatively  dry  to  wet.  After  achieving  dynamic  steady 
state,  the  y-ri  and  q-T)  trajectories  and  the  spectrum  of  the  response  time  series  are 
constructed.  The  shape-A  hillslope  and  boundary  condition  from  Figure  3.2  were  assumed 
for  the  simulations.  The  parameters  given  in  case  1  of  Table  3.3  were  implemented  in  the 
#  simulations.  Figure  3.20  illustrates  the  rainfall  forcing,  where  p(t)=fo-i-a  sin(2ic  t/T).  A 

diurnal  fluctuation  with  input  period  T  =  24  hours  was  assumed. 

Three  forcing  amplitudes  were  examined  as  described  in  Table  3.7.  The  amplitudes 
a  =  0.002  Ks  and  a  =  0.02  Ks  were  used  for  cases  IS  and  16,  respectively.  Case  17 
represents  a  =  0.  A  steady-state  solution  to  Richards  equation  based  on  the  rainfall 
®  intensity  of  0.001  Ks  was  the  assigned  initial  condition.  Note  that  all  three  cases  examined 

have  the  same  initial  condition  and  the  same  average  forcing  amplitude. 

Figure  3.21  shows  the  simulation  response  for  y,  t|,  and  q.  The  dotted  line  in  the 
figure  is  the  equilibrium  solution  based  on  fo  =  0.004  Ks  and  a  =  0. 


#  Figure  3.20. 


A  periodic  input. 


Table  5.4.  The  periodic  input  cases. 


case 

fc/Ks 

15 

periodic 

0.004 

0.002 

16 

ft 

ti 

0.02 

17 

II 

0.004 

0.0 

note:  Kg  =  0.013212  m/hr  =  saturated  hydraulic  conductivity  of  Guelph  Loam  soil. 
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It  was  observed  that  the  time  average  for  7  and  11  for  the  stationary  response  of 
cases  IS  and  16,  show  an  increasing  departure  from  the  equilibrium  case.  That  is.  as  the 
amplitude  of  fluctuation  increases,  the  mean  q  response  for  the  fluctuating  process  becomes 
increasingly  different  from  the  static  equilibrium  case.  Also  note  that  y,  X]  vary  inversely  as 
observed  earlier.  The  dynamic  equilibrium  for  the  two-state  variables  y  and  q,  depend  on 
the  amplitude  of  the  input  rather  than  the  mean  of  the  input,  demonstrates  a  fundamental 
nonlinearity  of  the  simulated  soil  water  process. 

Next,  we  look  at  the  output  behavior  at  dynamic  steady  state,  or  after  transient 
effects  have  died  out  Figures  3.22  and  3.23  show  the  time  series  and  their  corresponding 
power  spectra  for  case  IS  and  case  16,  respectively.  The  time  series  are  truncated  up  to  the 
lost  240  hours  of  the  signal.  The  power  spectra  were  calculated  based  on  a  record  length  of 
4096  hours  with  a  2  hour  sampling  interval  or  2048  sampling  points.  The  estimation  of 
power  spectrum  was  performed  using  MATHEMATICA  software.  The  unit  of  frequency 
is  expressed  in  terms  of  cycles/x,  where  1  is  the  record  length  of  4096  hours. 
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Figure  3.21.  Responses  for  y,  t],  and  q  for  cases  15-17. 
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Figure  3.22.  Time  and  frequency  domain  response  for  p,  y.  T),  r  and  q  for  case  15  (small  • 

amplitude  case  a  =  0.002Ks).  The  forcing  frequency  is  170.7  cycles/x,  t  =  4096  hours. 
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Figute  3.23.  Time  and  frequency  domain  responses  for  p,  y,  q,  r,  and  q  for  case  16  (large 
amplitude  case  a  =  0.02Ks).  The  forcing  frequency  is  170.7  cycles/t,  t  =  4096  hours. 
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spectral  analysis  of  the  time  series  is  used  to  determine  the  frequency  content  of  the 
time  series.  Recall  that  in  a  linear  system,  a  single  input  frequency  will  only  produce  the 
same  output  frequency.  On  the  other  hand  a  single  input  frequency  with  a  nonlinear  filter, 
produces  a  spectrum  of  frequencies  in  the  output.  Tne  fundamental  or  forcing  frequency  is 
the  largest  spectral  peak,  with  integer  multiples  of  the  forcing  frequency  indicating  the 
nonlinear  effect  This  form  of  nonlinearity  known  as  harmonic  distortion. 

Other  types  of  nonlinear  behavior  can  be  detected  via  signal  processing  tools. 
Higher  order  spectra  such  as  the  bispectrum  offers  the  possibility  of  detecting  the  role  of 
nonlinear  behavior  reflected  in  the  frequency  content  of  higher  time  moments  (e.g. 
skewness).  Our  initial  attempts  at  higher  order  spectral  analysis  for  stream-aquifer 
problems  in  this  research  effort,  are  descriiied  separately  in  Jin  and  Duffy  (1994). 

For  the  case  of  small  forcing  amplitude,  the  dominant  power  or  variance  in  y,  tj,  r 
and  q  spectra  is  concentrated  at  the  forcing  frequency  (Figure  3.22).  However  some 
variance  is  observed  at  harmonics  of  the  forcing  frequency.  The  peak  amplitude  in  the 
recharge  power  spectrum  is  3  times  as  large  as  the  peak  magnitude  in  rainfall  spectrum.  On 
the  other  hand,  the  peak  power  in  the  subsurface  outflow  spectrum  is  highly  damped 
compared  to  the  peak  amplitude  in  rainfall  or  recharge  power  spectra.  The  spectral  peak  for 
Y  and  q  spectra  are  nearly  the  same  magnitude. 

As  the  forcing  amplitude  increases  (case  16),  the  spectral  peaks  in  the  y,  tj,  r,  and  q 
series  also  increase  (Figure  3.23).  The  power  spectra  show  that  total  variance  in  the 
records  is  mainly  in  the  first  and  second  harmonics,  except  for  the  recharge  signal  where 
higher  harmonics  exist,  probably  related  to  the  estimation  problem  described  earlier. 
Clearly,  as  the  forcing  amplitude  increases,  the  amplitude  of  secondary  spectral  peaks  also 
increases.  As  in  case  15,  the  spectral  peaks  for  y  and  q  are  similar. 

Figure  3.24  shows  the  y-r\  phase  trajectories  constructed  from  the  stationary  time 
series.  It  is  found  that  the  phase  trajectories  cycle  in  a  clockwise  direction  but  the  looping 
is  not  very  broad.  Much  more  analysis  is  required  to  fully  understand  the  implications  of 
this  behavior. 

Figure  3.25  displays  the  q-q  phase  trajectories.  The  q-T|  phase  trajectories  change 
with  the  forcing  amplitude,  and  the  looping  structure  apparently  is  affected  by  the  secondly 
harmonics.  Again,  this  behavior  requires  additional  analysis  to  fully  understand  the 
structure. 
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Figure  3.23.  Phase  trajectories  for  for  cases  15  and  16. 
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Figure  3.24.  Trajectories  of  q-ii  for  cases  15  and  16.  The  direction  for  phase  trajectories 
is  clockwise  for  both  cases. 


3.5.4  Discussions  of  Results 

The  delayed  subsurface  outflow  peak  has  been  observed  through  simulation  with 
Richards  equation  under  a  specifled  hydrologic  conditions.  The  mathematical  origin  for 
delayed  subsurface  outflow  peak  appears  to  be  nonlinearity  inherent  in  soil  water  process. 
During  the  storm  event,  the  saturated  region  expands  upslope  due  rapid  infiltration  and 
recharge.  Following  the  storm  event,  accumulated  upslope  infiltration  moves  downslope 
primarily  as  a  saturated  layer  above  the  bedrock.  Since  this  lateral  subsurface  flow  takes 
some  time  to  reach  the  discharging  area  of  hillslope,  a  delayed  subsurface  flow  is  observed. 
This  mechanism  for  the  delayed  peak  has  little  to  do  with  the  expansion  of  the  saturation 
region  near  the  stream,  but  rather  depends  on  expansion  of  the  upslope  subsurface 
saturation  region.  Figure  3.25  illustrates  the  numerical  simulation  of  this  effect.  Each 
component  of  the  response  is  labeled.  The  magnitude  and  timing  of  the  delayed  subsurface 
outflow  peak  depends  on  factors  such  as  the  nature  of  input  forcing,  topography,  soil  type, 
and  storage  capacity  of  hillslope. 

The  dynamic  relationships  for  q(q)  based  on  Richards  equation  exhibit  hysteresis, 
even  though  the  soil  property  assumed  is  honiogeneous  and  non-hysteretic.  In  the  field. 
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Duffy  et  al.  (1991)  also  observed  a  hysteietic  relationship  between  baseflow  and  saturated 
water  storage  at  the  Upper  Sheep  Creek  hillslope  watershed.  They  also  provide  a 
qualitative  discussion  of  hysteresis  which  involves  the  interplay  between  slope  steepness, 
near-channel  storage,  the  distribution  of  inHltration  sources,  and  soil  heterogeneity. 
Without  any  complicating  factors  present  in  the  field,  hysteretic  behaviour  evident  in  the 
physics-based  model  under  the  idealized  conditions  appears  to  be  related  to  the  dynamic 
interaction  between  saturated  water  volume  and  subsurface  outflow. 


4.  FIELD  EXPERIMENT  AT  SHALE  HILLS 


4.1  Description  of  the  Experimental  Watershed  2 

The  Shale  Hills  Watershed  2,  a  7.9  hectare  forested  site  located  near  the  Penn  State 
campus  was  the  subject  of  an  NSF-funded  experimental  study,  initiated  in  the  1970's,  to 
evaluate  the  effects  of  antecedent  soil  moisture  on  stormflow  volumes  (J.  Lynch,  1976). 
The  site  is  a  forested  watershed  characterized  by  steep  slopes  and  narrow  ridges  typical  of 
the  low-lying  shale  hills  of  the  Ridge  and  Valley  Province  of  eastern  United  States.  The 
average  slope  gradients  are  16  and  18  degrees  on  the  south  and  north  aspects,  respectively. 
The  average  channel  gradient  is  4.5  %.  The  maximum  relief  is  61  meters  and  average  local 
relief  is  about  30  m. 

The  study  area  is  underlain  by  the  Rose  Hill  shale  about  213  m  thick  with 
occasional  interbedded  limestone.  Soil  depth  or  depth  to  bedrock  is  quite  uniform,  with  the 
ridges  only  slightly  thicker  than  the  valley  floor,  and  the  estimated  average  soil  depth  is  1.4 
meters  (Lynch,  1976).  Soil  types  on  the  watershed  are  Ashby  shaley-silt-loam  that  mantles 
the  moderate  to  steeply  sloping  valley  sides.  Blairton  silt-loam  is  found  on  the  gently 
sloping  upper  end  of  the  valley,  and  Ernest  silt  loam  underlies  the  lower  valley  floor. 
Leavesley  (1967)  and  Lynch  (1976)  provide  a  detailed  description  of  soil  physical 
properties,  the  soil  profile,  vegetation  characteristics,  and  history  of  land  use. 

At  Watershed  2,  Lynch  (1976)  instrumented  the  site  with  artificial  rainfall 
sprinklers,  and  simulated  rain  events  from  July  31, 1974  to  September  20, 1974.  During 
this  period,  8  sequential  artificial  storm  events  were  replicated  with  the  same  intensity  of 

6.1  mm  per  hour  and  duration  of  6.15  hours.  Rainfall  depth  was  measured  by  a  netwoik 
of  rain  gages:  six  recording,  40  standard  eight-inch  and  17  trough-type  nonrecording 
gages.  Stream  discharge  was  measured  at  4  different  weirs  with  continuous  recording 
water  level  recorders.  Soil  moisture  content  at  32  locations  was  monitored  using  neutron 
scattering  method.  At  each  site  soil  moisture  measurements  were  made  at  one  foot  intervals 
from  the  surface  to  bedrock.  Water  table  elevation  was  obtained  for  fifty  shallow 
observation  wells.  We  have  used  a  portion  of  this  data  base  to  investigate  a  method  for 
estimating  the  integrated  storage-flux  relationship,  as  part  of  a  longer  term  goal  of 
implementing  the  integral  balance  modeling  approach  at  the  site.  Table  4.1  below  gives  a 
summary  of  the  results  for  all  eight  storms.  Note  the  small  range  in  total  antecedent 
moisture,  and  the  order  of  magnitude  increase  in  total  storm  runoff. 
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antecedent 

total  storm 

rapid 

(telayed 

max.  peak 

Date 

soil  moisture 

runoff 

flow 

flow 

flow 

mm 

nun 

rrun 

mm 

m^s’l/km^ 

*8/01/74 

278.1 

4.55 

2.64 

1.90 

0.099 

*8/07/74 

287.5 

12.19 

8.36 

3.84 

0.201 

8/14/74 

290.1 

15.52 

11.94 

3.58 

0.253 

*8/19/74 

295.1 

15.90 

12.27 

3.63 

0.239 

8/23/74 

304.8 

20.75 

17.02 

3.73 

0.304 

8/27/74 

310.4 

23.34 

19.58 

3.76 

0.328 

*9/16/74 

313.2 

29.74 

25.78 

3.96 

0.526 

9/19/74 

316.2 

32.66 

28.88 

3.78 

0.628 

Table  4.1.  Stormflow  volumes  for  eight  6.1  mm/hr  artificial  storms  of  6.15  hour  duration. 
The  experiment  was  conducted  over  a  range  of  antecedent  moisture  conditions  (after 
Lynch,  et  al,  1976).  Note:  *  indicates  the  events  used  in  the  present  research. 


Figure  4.1.  illustrates  a  subset  of  the  Held  data  which  Lee  (1993)  has  analyzed  and 
these  results  are  presented  here.  Stream  stage  was  measured  continuously  during  and  after 
each  storm.  The  soil  moisture  and  groundwater  levels  were  measured  just  before  the  event, 
immediately  at  the  end  of  the  event,  and  daily  after  the  event,  until  the  runoff  response  was 
complete  or  a  new  rainfall  event  was  initiated.  Characteristic  curves  for  the  soils  are 
available.  Figure  4.1  also  shows  contour  lines  and  flow  paths  orthogonal  to  the  contours  at 
Watershed  2.  Six  soil  moisture  access  sites  and  eight  observation  well  sites  were  chosen 
for  this  study  between  weir  2  and  weir  4.  Four  events  were  selected  and  these  are  indicated 
in  Table  4.1.  Soil  moisture  and  water  table  elevation  measurements  were  available  before 
an  event  and  orrce  a  day  after  the  event 
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Figure  4.1.  Contour  and  path  lines,  and  monitoring  sites  at  Watershed  2. 

4.2  Defining  Similarity  Regions 

Duffy  et  al.  (1991)  have  proposed  a  local  affine  transformation  of  the  hillslope  as  a 
basis  scaling  state  variables  and  for  estimating  the  storage-flux  characteristics  of  a  hillslope. 
The  method  of  rescaling  the  trajectories  to  form  a  unit  hillslope  is  discussed  in  chapter  2. 
Of  course  the  critical  problem  is  to  determine  the  region  over  which  the  hillslope  geometry 
is  self-affine.  Figure  4.1  illustrates  a  hillslope  where  trajectories  were  found  to  be 
approximately  similar.  Recall  that  the  motivation  for  determining  similarity  regions  is  to 
rescale  the  land  surface,  bedrock,  water  table  and  soil  moisture,  such  that  sparse  and 
scattered  field  data  (water  levels,  bedrock  elevation,  etc.)  can  be  plotted  on  a  smooth  curve. 
Each  trajectory  starts  at  the  divide  and  truncates  at  the  stream  Qocal  base  level),  following 
the  steepest  path  orthogonal  to  contour  lines.  The  topography  scaled  in  this  manner  is 
locally  self-affine  since  two  length  scales  are  required,  one  for  elevation  and  one  for 
distance.  As  perfect  scaling  is  not  possible,  the  unit  hillslope  for  a  similarity  region  is 
constructed  in  a  statistical  sense,  by  fitting.  Once  the  similarity  region  is  determined  the 
hypsometric  relation  or  area  elevation  distribution  is  also  determined  so  that  spatial 
averaging  can  be  performed. 

By  trial  and  error  the  similarity  region  for  the  hillslope  in  Fig.  4.1  was  determined, 
and  the  fitted  surface  and  bedrock  topography  is  shown  in  Fig.  4.2.  The  difference 
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produces  the  ccmesponding  soil-colluvium  thickness  for  the  unit  hillslope.  The  fitted 
polynomials  along  with  the  rescaled  wator  table  at  a  particular  time  are  shown  in  Figure 
4.3,  which  together  define  the  saturated  fraction  of  the  unit  hillslope.  The  conditional 
hypsometric  distribution  for  the  hillslope  is  shown  in  Figure  4.4,  with  scaled  elevation 


O^rSl. 


Distance  Distance 

Figure  4.2.  Projection  of  the  surface  and  bedrock  trajectories  from  Figure  7  onto  the  unit 
hillslope.  Surface  shape  function:  Zo(s)=l  -  0.175  s  -  2.09  s^  +  1.27  Bedrock  shape 
function:  Z2(s)=  0.972  -  0.136  s-  2.411  s^  1.526  s^. 
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Figure  4.3.  Water  table  shape  function:  zi(s)=  0.829  +  0.374  s  -  2.807  s^  +  1.607  s^  and 
the  rescaled  unit  hillslope  zo(s)  with  corresponding  trajectory  of  bedrock  Z2(s).  Date:  19 
Aug.  1974. 


Figure  4.3.  The  area-elevation  relation  or  conditional  hypsometric  curve  for  the  hillslope  of 
Figure  4.1. 

In  the  field  experiment  (Lynch,  1976)  neutron  probes  were  used  to  estimate  the 
vertical  distribution  of  soil  moisture  at  1  foot  sampling  interval  down  to  the  bedrock  surface 
(e.g..  total  moisture  in  the  soil  column).  Water  table  elevation  measurements  were 
determined  independently  from  piezometers.  Integration  of  the  resulting  polynomials 
proceeded  as  outlined  in  chapter  2,  and  the  state  variables  for  spatially  integrated  soil 
moisture  and  saturated  storage  were  determined.  In  the  next  section  the  results  are  given 
for  four  simulated  rainfall  events  with  the  vaiying  antecedent  soil  moisture. 

4.3  Integrated  Storage-Flux  Response 

The  analysis  of  the  hillslope  study  for  the  four  events  are  shown  in  Figures  4.4-4.7 
Each  figure  shows  the  time  response  for  discharge  (Q  [m^/hr])  at  the  upstream  weir  4,  the 
downstream  weir  2,  estimated  recharge,  and  the  integrated  unsaturated  and  saturated 
storage  elements.  The  time  series  for  point  measured  soil  moisture  and  the  depth  to  water 
table  at  various  locations  on  the  hillslope  are  given  in  the  figure.  The  water  table  elevation 
and  depth-averaged  soil  moisture  are  plotted  versus  ground  surface  elevation. 
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Figure  4.4.  The  Q,  r,  y  and  T|  hydrographs,  and  time  and  space  distribution  for  depth- 
averaged  moisture  and  water  table  for  8/01/74  storm.  The  ground  surface  and  water  table 
refer  to  elevation.  The  position  for  moisture  access  tube  and  observation  well  locations  is 
shown  in  Figure  4.1. 
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Figure  4.5.  The  Q,  r,  y  and  T|  hydrographs,  and  time  and  space  distribution  for  depth- 
averaged  moisture  and  water  table  for  8/7/74  storm.  The  ground  surface  and  water  table 
refer  to  elevation.  The  position  for  moisture  access  tube  and  observation  well  locations  is 
shown  in  Figure  4.1. 
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Figure  4.6.  The  Q,  r,  y  and  t)  hydrographs,  and  time  and  space  distribution  for  depth- 
averaged  moisture  and  water  table  for  %l\9nA  storm.  The  ground  surface  and  water  table 
refer  to  elevation.  The  position  for  moisture  access  tube  and  observation  well  locations  is 
shown  in  Figure  4.1. 
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Figure  4.7.  The  Q,  7  and  hydrographs,  and  time  and  space  distribution  for  depth- 
averaged  moisture  and  water  table  for  9/16/74  storm.  The  ground  surface  and  water  table 
refer  to  elevation.  The  position  for  moisture  access  tube  and  observation  well  locations  is 
shown  in  Figure  4.1. 
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The  discharge  hydrographs  at  downstream  weir  2  display  both  single  and  double 
peaks,  while  no  delayed  streamflow  peak  is  observed  at  upstream  weir  4.  The  double 
peaks  are  present  in  the  8/01/74  and  8/19/74  storms,  a  minor  second  peak  occurs  on 
8/7/74,  and  the  9/16/94  storm  has  no  double  peak.  For  8/01/74  storm,  the  delayed  peak 
discharge  is  smaller  than  the  initial  peak,  but  the  second  peak  is  larger  than  the  first  for 
8/19/74  storm.  The  presence  of  double  peaks  spears  to  vary  with  antecedent  soil  moisture 
conditions,  and  wetter  initial  conditions  cause  the  second  peak  to  diminish  or  disappear. 
We  note  that  the  discharge  at  weir  4  is  small  in  magnitude  and  lags  discharge  at  weir  2. 

Recharge  hydrographs  are  constructed  based  on  an  explicit  finite  difference  form  of 
the  integral  balance  equation  for  the  saturated  zone  as  discussed  in  section  3.5,  where 
discrete  q  and  r\  time  series  are  used.  During  the  event,  recharge  rapidly  increases,  while 
post-event  recharge  decreases  with  a  sign  change  from  positive  to  negative  during  the 
drainage  period.  The  r  is  positively  correlated  with  q  and  negatively  correlated  with  y 
except  for  8/1/74  storm  (dry  conditions).  Saturated  storage  increases  for  each  storm  event, 
and  decreases  during  the  recession  period.  Unsaturated  storage  y  shows  an  inverse  relation 
to  q  response.  The  distribution  of  soil  moisture  and  water  table  elevation  measured  at  3 
locations  demonstrate  the  difference  between  the  point  measured  response  and  the 
integrated  response.  Note  that  at  a  point  on  the  hillslope  the  soil  moisture  response 
increases  during  and  then  decreases  after  the  storm,  while  the  integrated  soil  moisture 
decreases  and  then  increases.  It  is  important  to  emphasize  that  the  competitive  nature  of  y 
and  q  dynamics  observed  in  the  numerical  experiments  is  also  observed  in  the  field  data. 
Another  interesting  observation  is  that  the  overall  range  of  integrated  soil  moisture  change 
is  quite  small,  less  than  4  %  by  volume  for  the  entire  experiment. 

Although  the  data  are  sparse  for  all  storms  examined,  the  time  scales  for  discharge 
and  storage  are  different  when  comparing  the  storm  event  with  the  drainage  period.  If  we 
deHne  the  time  scale  as  the  slope  of  the  curve  Q  versus  dQ/dt,  then  during  the  storm  phase 
the  time  scale  is  seen  to  be  much  shorter  than  during  recession  phase.  This  may  also  be  the 
result  of  the  competitive  relationship  between  capillary  and  gravity  storage  during  the  storm 
and  recession  periods. 

The  relationship  between  altitude  and  water  table  elevation  is  evaluated  at  two  times, 
before  and  after  the  storm.  The  relationship  shows  that  the  water  table  is  highly  correlated 
with  surface  topography.  This  empirical  correlation  provides  support  for  the  unit  hillslope 
scaling  method  for  water  balance  studies.  The  depth-averaged  soil  moisture  is  high  near  the 
'stream  as  expected,  is  lower  over  the  mid-slope  region,  and  interestingly  is  somewhat 
higher  again  in  the  ridge-top  region.  Again,  a  clear  pattern  exists  within  the  similarity 
region. 
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Figure  4.8  shows  the  and  q-ii  trajectories  for  all  storms  analyzed  in  this  chapter. 

The  trajectories  exhibit  an  inverse  relationship  which  indicates  a  competitive  relation 

between  the  state  variables.  For  the  dry  initial-moisture  condition  of  the  8/01/74  event,  y  # 

increases  with  Ti  dining  the  storm,  while  y  increases  and  T)  decreases  during  the  recession 

phase.  The  looping  behavior  is  seen  the  y-r;  phase  portrait,  but  the  width  of  the  loop  is 

small.  The  looping  or  hysteretic  behavior  is  much  more  pronounced  in  q-r^  trajectories. 

The  q-Tl  paths  reveal  clockwise  evolution  for  all  storm  events  and  its  shape  and  magnitude 
vaiy  with  antecedent  soil  moisture  condition.  . 

Figure  4.9  illusuates  the  y(ii)  and  q(ii)  relationship  for  the  data  of  all  four  events.  A 
linear  function  was  fitted  to  the  relationship  y(Ti).  and  a  quadratic  polynomial  was  fit  to  the 
relationship  q(Ti). 

Figure  4. 10  illuiitrates  the  structure  of  r-y-T)  trajectories  for  each  storm.  The  figure  • 

also  indicates  the  projection  of  r-y-t]  phase  trajectories  onto  orthogonal  planes  of  r-y,  r-Tj, 
and  Y-T|.  Although  the  data  are  sparse,  recharge  is  negatively  correlated  with  unsaturated 
storage  for  the  recessions  of  the  8/7/74  and  8/19/74  events,  and  the  trajectory  in  the  r-y 
plane  shows  a  counter  clockwise  motion. 
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Figure  4.8.  The  y-T|  and  q-T)  phase  paths.  Numbers  denote  the  position  of  phase  trajectory 
in  hours  from  the  start  of  rainfall  simulation. 
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4.4.  Discussion  of  Results 

The  single  most  important  result  in  this  analysis  of  the  Shale  Hills  data,  is  the 
observed  inverse  relationship  for  the  phase-plane  trajectory  of  unsaturated  and  saturated 
storage,  y-ii.  This  inverse  relationship,  also  found  in  the  numerical  experiments,  is  the 
simple  result  of  competition  between  y  and  t|  for  storage  volume  on  the  hillslope.  The 
forces  governing  storage  in  the  unsaturated  state  (capillarity  and  gravity),  and  the  forces 
governing  saturated  stora(*e  (gravity),  introduce  different  time  scales  of  response  in  each 
zone. 

The  dynamic  constitutive  relationship  q(ti)  manifests  hysteretic  (looping)  behavior 
with  a  clockwise  rotation  during  storm  period.  A  quadratic  polynomial  for  q(r|)  was 
determined  by  regression  analysis,  and  this  model  should  be  useful  for  modeling 
nonhysteretic  subsurface  outflow  and  saturated  water  storage.  However,  the  role  of 
hysteresis  appears  to  be  significant. 

It  is  shown  that  local  similarity  of  the  hillslope  shape  provides  a  useful  scaling  of 
scattered  geomorphic  and  hydrologic  data.  Within  the  conditional  similarity  region  the 
water  table  surface  is  projected  onto  the  unit  hillslope  geometry,  as  is  the  soil  moisture. 
This  projection  allows  the  rescaled  water  table  and  soil  moisture  to  be  fit  to  a  smooth  curve, 
which  is  integrated  over  the  hiUslope. 

The  different  time-scales  for  rising  and  recession  stages  offer  some  insight  into  the 
problem  of  sampling  frequency  in  the  field  experiment.  The  quick  response  during  the 
storm  event  requires  greater  sampling  resolution  than  the  drainage  period.  Comparison  of 
the  response  time-scales  of  storage-flux  from  the  Richards  equation  and  from  the  field  data, 
show  that  the  watershed  responds  much  more  rapidly  than  the  numerical  model.  Since  the 
soil  types  used  in  the  numerical  study  are  not  much  different  from  the  soils  encountered  in 
the  field,  the  large  difference  in  characteristic  times  may  very  well  be  due  to  heterogeneity. 
The  effective  hydraulic  conductivity  of  the  watershed  is  much  higher  than  would  be 
determined  by  the  texture  of  the  individual  soil  types.  This  would  suggest  that  the 
integrated  performance  of  the  hillslope  or  watershed,  and  the  'effective'  soil  hydraulic 
properties  cannot  be  determined  by  laboratory  (small)  samples  but  must  be  determined 
directly  in  the  field.  We  propose  to  investigate  this  in  future  work. 
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5.  THE  TWO-STATE  MODEL 


5.  THE  TWO-STATE  MODEL 

5.1  The  Two-State  Model  and  its  Phase  Space 

The  numerical  experiments  with  the  Richards  equation,  and  the  inteipretation  of  the 
Shale  Hills  field  experiment,  offer  reasonable  confirmation  that  spatially  integrated  state 
variables  can  describe  the  dynamics  of  the  rainfall-recharge-runoff  process.  In  this  section 
an  initial  attempt  is  made  at  formulating  the  particular  dynamical  system.  Although  the 
numerical  experiments  and  the  field  data  suggest  that  nonlinear  constitutive  relations  are 
ultimately  necessary  to  represent  subsurface  dynamics,  it  is  always  useful  to  begin  with  the 
simplest  linearized  form.  For  a  linear  two-state  model,  solutions  are  found  for  competitive 
and  noncompetitive  state  variables,  responding  to  varying  initial  conditions  and  external 
forcing  (rainfall).  Recall  that  the  two-component  state-space  form  of  the  model  is  given  by 

^=i(p,e.t)-r(n,Y,t)  (41) 

dt 

^=r(Ti,7,t)-q(Ti,t)  (42) 

dt 

where  i(p,e,t)  is  surface  infiltration  rate,  or  that  flux  which  infiltrates  the  soil  less 
evaporation,  r(Tj,y,t)  is  the  recharge  rate  per  unit  surface  area  to/from  the  water  table,  and 
q(T)>t)  the  volumetric  subsurface  runoff  per  unit  area. 

The  dynamic  behavior  of  the  two-state  system  is  graphically  represented  by  the 
state-space  or  time  series  of  7(t)  and  Ti(t),  and  by  constructing  the  trajectory  of  a  parametric 
curve  Y(t)-Ti(t)  referred  to  as  the  phase  space  or  phase  portrait.  Figure  5.1  illustrates  the 
two  represemations  of  the  state  variables.  The  parametric  curve  y-q  defines  the  motion  or 
trajectory  of  the  state  variables  with  time  implicit  for  a  given  initial  condition  (y^,  q^),  and 

is  referred  to  as  the  phase  portrait.  All  possible  trajectories  would  define  the  phase  space. 
Both  have  been  used  in  this  report.  The  next  step  is  to  formulate  the  specific  constitutive 
relations. 
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Figure  5.1.  Representation  in  terms  of  state-space  and  phase-space. 


5.2  Linear  Constitutive  Relations  and  State  Space 


A  general  two-state  linear  model  has  the  form 
dX] 
dt 
dx., 


—  u,  +  a,[Xi  UjjXj 


dt 


=  Uj+a2iX, +3^2X2 


(44) 


The  input  vector  will  be  simplified  in  the  present  case{ui,0},  since  rainfall  ui  is  the  only 
input  to  the  system.  First  the  state  variables  must  also  be  modified  to  include  residual 
storage,  xi=(y-^o)  and  X2=(ti-i1o).  Residual  storage  in  the  case  of  the  unsaturated  zone  is 
that  depth  of  water  which  is  held  in  storage  by  capillary  forces  after  gravity  drainage  is 
complete.  The  residual  storage  for  the  saturated  region  is  that  depth  of  water  which  exists 
below  local  base  level.  The  task  of  determining  the  coefficients  aij  for  the  two-state  model, 
relies  on  determining  a  series  of  rate  constants,  the  form  of  which  are  determined  from  the 
numerical  and  field  experiments. 


5.3  Noncompetitive  Model  for  y-Ti 


The  first  parameterization  makes  the  assumption  that  outflow  from  each  storage 
reservoir  is  proportional  to  the  remaining  storage  in  that  reservoir  and  that  transfer  of 
moisture  from  one  state  to  the  other  is  also  directly  proportional.  The  latter  condition 
implies  that  increasing  moisture  in  the  unsaturated  zone  should  increase  recharge  to  the 
saturated  region.  Since  there  is  only  a  single  output  to  this  linear  system  (neglect 
evaporation  for  the  time  being),  the  output  (subsurface  flow  to  the  stream)  is  linearly 
proportional  to  the  storage  in  the  saturated  reservoir.  In  standard  notation,  the  coefficients 
are  defined  in  terms  of  the  rate  constants  kij 


(45) 


a,,  -  -kj, 
^21  “  ^21 


^12  ~  ^12 


^22  “  ~(^02  ^12  ) 

where  the  subscript  on  the  rate  constants  {i.j}  represent  the  direction  and  origin  of  the 
transfer  of  fluid  respectively.  The  output  rate  constant  from  the  second  reservoir  is  ko2. 
So  this  form  of  the  model  takes  the  form 


^  =  P~k2,(y-Yj+k,j(Ti-TiJ 


dt 


(46) 


^2,  (Y  -  Yo )- (ko2+ 1^12  )(T1- Tie) 


We  can  refer  to  this  formulation  as  'noncompetitive',  since  increases  in  one  state  variable 
serves  to  enhance  the  other.  The  equilibrium  or  steady-state  is  found  by  setting  the  time 
derivative  to  zero 


Y  =  Y„+i^  +  ^(Tl-tlo) 

kj,  kj, 

k  -Hk 

Y  =  Yo-*-*^\  (Ti-n,) 

Kj] 

The  direct  proportionality  of  y-r\  clearly  shows  the  noncompetitive  nature  of  this 
pa’  ameterization  of  the  model.  The  single  equilibrium  or  fixed  point  is  given  by 


fk  +k  ' 

'y(-)=p 

V  *■21*^02  > 


+  Yo 


(48) 


Tl(oo)  = 


+  T1<, 


The  general  form  of  the  solution  for  (46)  v'hen  p=^  is  given  by 

x(t)  =  a.e*^‘‘ci  +  a^e^^'icj  ( 53 ) 

where  (xi.xi)  are  the  state  variables,  {ai.ai}  are  scalars  determined  from  initial 
conditions,  {CpCj}  are  the  eigenvectors  of  the  coefficient  matrix  (45),  and  {Xi,  X2}  are  the 
corresponding  eigenvalues.  In  this  case  we  require  that  be  real  and  negative  for  the 
solution  to  be  stable. 

The  phase  portrait  for  a  range  of  initial  conditions  and  the  indicated  parameters  is 
shown  in  Figure  5.2.  The  phase  portrait  represents  all  possible  solutions,  including  some 
that  are  clearly  not  physically  plausible.  The  time  series  for  a  particular  initial  state  is 
shown  in  Figure  5.3.  In  the  next  section,  the  results  of  the  noncompetitive  model  will  be 
compared  with  the  case  where  soil  moisture  and  saturated  storage  compete  for  the  same 
hillslope  storage  volume. 
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Figure  5.2  The  phase  portrait  for  the  noncompetitive  model  with  the  following  parameters; 
ki2=0.1,  k2i=0.1,  ko2=0.2,  p=0.02,  Yo=0-7  and  110=0.1. 


Figure  5.3  The  time  series  for  same  parameters  as  in  Figure  5.2  but  with  initial  conditions 
{Yo=.6,i1o=.4}. 

#  5.4  Cor  "'(’titi  e  Model  for  y-Ti 

Introducing  competition  between  state  variables  simply  means  that  increasing  the 
storage  in  one  reservoir  will  cause  a  decrease  of  storage  in  the  second  reservoir.  Recall  that 
^  both  the  numerical  experiments  and  the  Shale  Hills  field  data  exhibit  an  inverse  relationship 

for  y-tI.  except  during  very  dry  conditions.  The  concept  of  competition  may  be 
counterintuitive  at  first,  until  we  recall  that  the  state  variables  represent  the  average  or  total 
moisture  in  each  reservoir.  Under  well-drained  conditions,  the  saturated  storage  is  at  its 
lowest  point,  t|->t|o,  while  the  moisture  store  may  be  at  a  high  level.  During  a  storm, 
^  infiltration  increases  the  soil  moisture  at  first,  but  at  some  later  time,  the  water  table  begins 

to  capture  the  soil  moisture  volume.  Below  the  root  zone,  the  initial  unsaturated  volume 
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contains  substantial  moisture.  Thus,  a  small  amount  of  recharge  may  produce  a  rapid 
increase  in  the  water  table,  and  effectively  reduce  the  volume  of  moisture  in  storage  in  the 
unsaturated  zone.  However,  both  the  numerical  and  field  experiments  suggest  that,  during 
the  early  part  of  the  event,  both  reservoirs  increase.  The  transition  from  noncompetitive  to 
competitive  state  variables  is  an  important  concept,  which  we  examine  next 
The  coefficients  are  again  defined  in  terms  of  rate  constants  ky 

( 49 ) 

*21  ”  ^21  *22  ”  “(^02  ”  ^12  ) 

Note  that  this  version  of  the  model  is  only  different  in  that  k\2  has  changed  sign.  The 
dynamical  model  for  competition  between  state  variables  takes  the  form 


dt 


dt 


(50) 


k2,(Y-Yo)-(ko2-k,j)(Ti--n„) 


The  equilibrium  or  steady-state  is  again  found  by  setting  the  time  derivative  to  zero 

P  k,, ,  . 

Y  =  Yo+T^-7^(t1-t1„) 


Y  =  Yo  + 


*■21  *■21 
l^-k 


(51) 


^(T1-T1„) 


‘•21 


The  single  equilibrium  or  fixed  point  for  this  linetu  system  is  given  by 


Y(«)  =  p 


^02  ~  ^12 
kjjkflj 


+  Yo 


(52) 


T1(oo)  =  -^  +  T1o 


The  general  solution  for  (50)  when  p=0  is  given  by 

x(t)  =  e®'  cos(qt)  +  Cj  sin(qt)}  ( 53  )  • 

where  the  eigenvalues  in  this  case  may  be  complex  Xi=o±iq,  and  the  coefficients  (Ci  .C2 ) 
are  related  to  the  eigenvectors  of  the  coefficient  matrix  (49).  As  long  as  the  real  part  of  Xj  is 
negative,  the  solution  is  stable.  The  reader  is  referred  to  Beltrami  (1987)  for  details  of  the 
dynamic  modeling  approach.  ^ 

The  phase  portrait  for  a  range  of  initial  conditions  and  the  indicated  parameters  is 
shown  in  Figure  5.4,  and  we  see  that  the  solution  takes  the  form  of  a  spiral  as  it  approaches 
the  single  equilibrium  point.  The  phase  portrait  represents  all  possible  solutions  for  the 
given  parameter  set,  including  some  that  are  not  physically  plausible.  It  should  be  noted, 
that  the  addition  of  nonlinear  constitutive  relations  will  be  necessary  to  constrain  the  ^ 

nonphysical  trajectories,  or  those  where  the  state  variable  is  negative.  The  time  series  for  a 
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particular  initial  state  of  the  competitive  model  is  shown  in  Figure  S.S.  Note  the  substantial 
change  in  response. 

Although  our  numerical  and  field  experiments  show  that  nonlinearity  is  the  rule 
rather  than  the  exception,  the  linear  models  outlined  here  provide  a  first  step  towards 
understanding  the  physical  mechanisms  operating  as  they  are  expressed  by  integral 
variables  of  the  hillslope  system.  Future  work  will  explore  in  much  greater  detail  the 
general  nonlinear  form  of  the  dynamical  model. 


Figure  5.4  TTie  phase  portrait  and  equilibrium  or  fixed  point  for  the  competitive  model  with 
the  following  parameters:  ki2=0.1,  k2i=0.1,  ko2=0.2,  p=0.02,  Yo=0.7  and  rio=0.1. 


Figure  5.5  The  time  series  for  the  competitive  model  with  same  parameters  as  in  Figure 
5.4  but  with  initial  conditions  {7o=-6,'no=.4}. 
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6.  DIGITAL  ELEVATION  DATA  AND  TOPO-HYDROLOGIC 

STRUCTURES 


6.1  Fractal  Topography  and  Similarity  Regions 


As  part  of  this  research  an  initial  attempt  was  made  to  examine  similarity  of  digital 
topography,  and  the  scale  over  which  similarity  might  be  expected.  It  has  been  proposed 
that  scale  invariance  of  random  topographic  surfaces  represents  a  fundamental  characteristic 
of  landform  texture  and  thus  classification  (Turcotte,  1992,  Mandelbrot,  1985,  others). 
Digital  topography  for  two  7.S  min.  quadrangles,  one  in  the  Colorado  Plateau,  and  one  in 
the  Appalachian  Plateau  were  used.  As  part  of  the  analysis,  algorithms  were  developed  to 
analyze  the  data,  and  these  are  described  below.  Although  the  analysis  was  limited  in 
scope,  the  results  point  to  a  critical  avenue  for  future  research. 

An  elementary  characterization  of  random  surface  texture  requires  two  parameters, 
the  relative  roughness  or  variance  of  the  surface,  and  the  distribution  of  variance  with  map 
scale.  For  a  self-affine  surface  z(x,y),  these  properties  are  measured  by  the  root-mean- 
square  roughness,  Cz  and  the  exponent  P  of  the  variance  spectrum  r(k;P),  where  the 
spectrum  has  the  form 

r(k)  -  k-P  ( 54 ) 

For  this  model  of  topography,  the  parameter  P  represents  the  way  that  variance  is 
distributed  with  wave  number  k  or  wavelength,  k'i=X.  Oz  and  P  may  also  be  thought  of  as 
the  scale  and  shape  parameters  respectively,  of  the  correlation  structure  in  z(x,y).  A  useful 
property  of  the  spectrum  is  that  integration  over  wave  number  k  defines  the  overall  variance 
of  the  surface,  or  in  terms  of  the  rms  roughness 


where  ko  is  the  smallest  and  ki  is  the  largest  wave  number  to  exhibit  self-affine  scaling.  A 
relation  between  roughnes  Oz  and  map  length  can  be  determined  by  integrating  (54)  (Power 


and  TulUs,  1988) 


r  o  V'* 

*  U-iJ  “ 


(56) 


where  a  is  a  constant,  and  Lo=:l/ko  is  the  longest  wavelength  over  which  the  scaling 
behavior  is  observed.  The  fractal  dimension  D  for  the  surface  can  be  expressed  in  terms  of 
the  spectral  slope  P 

D  =  I^.  (57) 


74 


A  spectral  estimation  algorithm  was  developed  in  MATLAB  based  on  the  2D  fast 
fourier  transform.  The  field  was  assumed  to  be  statistically  isotropic,  and  the  radial- 
symmetric  spectral  estimate  was  determined.  By  plotting  Inf  versus  Ink  for  a  subset  of 
the  DEM,  the  spectral  slope  and  fractal  dimension  (57)  were  estimated  and  a  comparison  of 
the  distribution  of  values  were  made.  In  this  example  two  7.5  minute  USGS  DEM's  are 
examined,  one  located  in  the  Colorado  Plateau  near  Price,  Utah  and  a  second  in  the 
Appalachian  plateau  in  northwestern  Pennsylvania.  Each  spectrum  estimate  was 
determined  for  a  120x120  grid  subset  of  the  DEM,  and  the  results  are  shown  in  Figure  6.1. 
The  approximate  upper  and  lower  cutoff  frequencies  are  indicated  and  the  mean  fractal 
dimension  D  and  ^  is  given.  As  the  grid  interval  is  approximately  30  meters  the  maximum 
length  scale  of  similarity  k:‘is  Lo~2km,  and  the  minimum  length  scale  kj^'  appears  to  be 
on  the  order  of  Li-lOOm. 

Although  this  is  only  a  preliminary  analysis  based  on  a  small  amount  of  digital 
topography,  the  fractal  dimension  and  spectral  slope  really  do  not  vary  all  that  much  for 
these  cases.  Future  work  will  expand  this  phase  of  the  study  and  in  particular  to  examine  5 
m  digital  data.  Understanding  the  transition  of  similarity  relationships  from  the  hillslope  to 
watershed  scale  will  require  greater  detail  in  the  topographic  data.  A  future  goal  is  to 
examine  multi-scale  hypsometric  relationships  from  the  digital  topography. 
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Figure  6.1.  The  estimated  spectra  for  Slate  Run  PA.  and  Wellington,  UT 
for  7.5  minute  quadrangle.  Each  trace  represents  a  120x120  subgrid  from  the 
DEM. 


6.2  Application  to  Climate  Variables  in  Mountainous  Regions 

In  chapter  2  we  considered  a  general  space-time  representation  for  the  state  variable 
y(X,t}  ofthefoim 

y(x.t)  =  ^yk(t)ak(*)  ( 3 ) 

k-l 

where  yk(t)  represents  the  kn,  amplitude  and  a^(k)  are  the  corresponding  spatial  shape 
functions  with  respect  to  the  vector  of  independent  variables  or  coordinates  x={  xi,X2>X3 } . 
From  (3)  the  problem  is  to  determine  the  ‘effective'  dimension  of  the  independent 
variables,  and  their  corresponding  shape  functions.  A  good  example  is  the  distribution  of 
precipitation  and  temperature  in  mountainous  regions,  where  the  climatic  variables  are  often 
reduced  to  simple  functions  of  elevation.  Let's  say  that  our  state  variable  y  is  a  function  of 
two  independent  coordinates,  elevation  and  aspect  or  azimuth  and  further  that  y  is  a 
stationary  process  in  dme.  The  process  y(<t>,8,t)  is  then  decomposed  into  a  dme  average 
(y)  and  a  stationary  time-fluctuation  (y') 

y(<l>,e,t)=y«j),e)+y'((|),e,t)  (58) 

where  the  time  average  y((t>,6)  is  a  function  of  elevation  (<|>)  and  aspect  (6),  while  the 
perturbation  y*  (<(>,6,t)  varies  with  elevation,  aspect  and  time.  Now  let  y',  have  the 
expansion 

y'(<(>,t)  =  X  (59) 

I: 

where  y^Ct)  is  a  zero  mean  stationary  process  in  time  representing  stationary  harmonic 
and/or  random  contributions  to  the  overall  process,  and  represents  the  basis 

functions  or  shape  functions  for  the  same  component  of  the  y  fluctuation.  The  y-process 
now  has  the  representation 

y«t».e,t)  =  y«|»,e)+Xyk(t)a,((j>,e)  .  (60) 

k 

The  hydrologic  problem  is:  How  does  one  determine  the  shape  functions  ak((|>,6)  from 
time  series  measurements  of  y? 

As  an  application  of  (60),  Fan  and  Duffy  (1993)  applied  the  decomposition  (60)  the 
seasonal  structure  of  monthly  temperature  (T)  and  precipitation  (P)  records  for  the  Wasatch 
Front  area  of  north-central  Utah.  Each  time  series  was  expressed  as  the  sum  of  a  long  term 
mean,  a  seasonal  cycle,  and  a  residual  stochastic  process.  In  this  example  the  shape 
functions  were  assumed  to  be  only  functions  of  elevation,  as  all  climate  stations  were 
located  on  the  western  slopes  of  the  Wasatch.  They  then  examined  the  elevation 
dependence  of  the  mean,  the  amplitude  and  phase  of  Fourier  harmonics  associated  with  the 
seasonal  cycle,  and  the  variance,  autocorrelation  and  skew  of  the  residual  process.  Based 
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on  coirelation  with  elevation  and  each  of  the  above  parameters,  a  simple  parametric  model 
of  T  and  P  as  a  function  of  elevation  (z)  and  season  (t)  is  developed,  which  preserves  the 
seasonal  variability  in  the  data.  The  model  provided  an  empirical  method  for  estimating 
spatially-distributed  and  time-variable  climatic  input  patterns  from  point  observations.  The 
application  demonstrated  that  seasonal  trajectory  of  the  P(z,t)-T(z,t)  phase  portrait  can  be 
diagnostic  of  synoptic  climate  conditions  for  interpretation  of  large  seasonal  variations 
about  the  normal  mean  conditions.  Finally,  the  model  is  applied  to  digital  elevation  data 
(DEM)  of  the  Wasatch  Front  area,  to  generate  seasonal  animations  of  temperature  and 
precipitation  fields.  The  emphasis  in  this  study  is  placed  on  the  space-time  structure  in 
historic  records  and  proposing  an  empirical  method  to  incorporate  such  structures. 

6.3  Estimating  Precipitation-Temperature  Shape  Functions 

One  expects  the  mean  temperature,  px.  to  decrease  with  elevation  in  a  manner  close 
to  the  lapse  rate,  and  the  mean  precipitation,  pp,  to  increase  with  elevation  as  a  result  of  the 
orographic  effect.  Figure  6.2  and  6.3  illustrate  monthly  time  for  temperature  and 
precipitation  at  three  elevations,  along  with  their  spectra.  Fan  and  Duffy  attempted  to 
determine  whether  higher  statistical  moments,  namely  tlie  variance  and  skew,  might  also  be 
functions  of  altitude.  The  variance  of  the  climate  time  series  was  further  decomposed  into 
the  contributions  from  seasonal  harmonics.  The  state  variables  have  the  form 

m 

T(z,t)  =  PT(z)  +  ^aT,(z)  cos{2m/Jt  +  <|>T,(z)])  +  N[OT(z).PT(zX^T(z)]  (61 ) 

id 

m 

P(z,t)  =  Pp(z)  +  ^ap^(z)  cos{27ti/Jt  +  (pp^(z)]}  +  N[ap(z).pp(z),Op(z)]  (62 ) 

i>l 

where  the  noise  parameter  is  a  function  of  the  residual  variance,  the  spectral  slope,  and  the 
skew  coefficient  respectively.  Each  parameter  in  (61)  is  estimated  from  the  data  and  plotted 
against  elevation  as  shown  in  Figure  6.4.  Note  that  each  data  point  represents  a  climate 
record.  The  upper-left  plot  indicates  that  monthly  mean  temperature  decreases  with 
elevation  at  a  rate  of  6.06°C  /  1000m.,  slightly  lower  than  the  sea-level-average  lapse  rate 
of  6.5oC  /  1000m.  The  upper-right  plot  gives  the  amplitude  of  the  first  three  harmonics 
used  to  describe  the  seasonal  cycle.  Three  harmonics  are  found  to  be  the  maximum  number 
needed  to  meet  the  criteria.  The  plot  shows  no  elevation  dependence  for  the  second 
harmonic  (f=l/6  cpm,  cycle  per  month).  For  the  first  (f=l/12  cpm)  and  the  third  (f=l/4 
cpm),  the  amplitude  decreases  with  elevation,  or  at  higher  elevations  one  observes  less 
seasonal  difference,  and  the  air  stays  relatively  cool  all  year  round.  The  larger  seasonal 
amplitude  at  low  elevation  may  also  reflect  the  winter  inversion  which  causes  cooler  winter 
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temperature.  The  middle-left  plot  shows  a  lack  of  elevation  dependence  for  all  three  phase 
components.  On  the  average,  the  peak  occurs  in  the  ICi  month  (July,  warmest  month)  for 
the  first  harmonic,  the  5th  month  (May)  for  the  second,  and  between  the  3rd  and  4th  month 
Gate  March  to  early  April)  for  the  third.  The  variance  of  residual  noise  is  given  in  the 
middle-right  plot,  and  it  shows  large  scatter  and  a  weak  correlation  with  elevation.  The 
lower-left  plot  gives  the  distribution  of  the  absolute  values  of  spectral  slope  for  the  residual 
noise.  The  lower-right  plot  shows  the  elevation  distribution  of  the  skew  coefficient  for  the 
noise  process,  which  indicates  cooler  (than  the  mean  T)  in  the  valley  and  warmer  (than  the 
meanT)  at  higher  altitudes. 

Figure  6.5  illustrates  the  elevation  distribution  of  precipitation  moments.  The 
upper-left  plot  is  the  long  term  mean  monthly  precipitation,  fitted  to  a  linear  variation  with 
elevation.  The  upper-right  plot  gives  the  amplitude  of  the  first  five  harmonics  used  to 
describe  seasonal  oscillations.  The  residual  contribution  to  the  variance  is  higher  for 
precipitation  than  for  temperature  (0.1%).  Five  harmonics  are  found  to  be  the  maximum 
necessary.  The  plot  shows  th^*  only  the  first  (f=l/12  cpm)  and  the  third  (f=l/4  cpm) 
harmonics  have  a  significant  elevation  effect.  The  general  increase  of  seasonal  fluctuation 
with  altitude  can  be  explained  as  follows  (Williams  and  Peck,  1962).  In  the  Great  Basin 
region,  precipitation  in  winter  is  mostly  of  frontal  type,  resulting  in  a  strong  orographic 
signature.  In  the  summer  months,  precipitation  is  primarily  from  local  convective 
thunderstorms,  without  apparent  elevation  effects.  At  the  watershed  scale,  this  leads  to  a 
small  seasonal  difference  in  precipitation  at  low  elevation,  and  a  large  seasonal  difference  at 
high  elevation.  The  middle-left  plot  shows  the  phase  relationship  for  the  five  harmonics. 
No  significant  altitude  correlation  is  found  except  for  the  first  harmonic  (f=l/12  cpm), 
whose  peak  comes  between  the  12th  and  the  1st  month  (December  to  January),  the  season 
when  frontal  storms  from  the  North  Pacific  are  most  active.  The  second  harmonic  (f=l/6 
cpm)  appears  to  be  related  to  a  large  scale  low-pressure  cell  called  the  Utah-Nevada  Low 
(Jeppson,  et  al.,  1968)  that  occurs  in  the  late  spring  months  (April  to  May)  and  the  mid-fall 
season  (October),  a  period  of  widespread  thunder-storms.  The  middle-right  and  lower-left 
plots  give  the  variance  and  the  spectral  slope  b  of  the  residual  series.  It  is  observed  that  the 
variance  and  the  persistence  of  the  residual  stochastic  process  increase  with  altitude.  The 
lower-right  plot  shows  the  large,  positive  skew  coefficient  in  the  residual  noise  at  all  sites 
and  the  lack  of  elevation  dependence.  Large  positive  skew  coefficient  is  a  characteri.<;tic  of 
many  hydrologic  processes  in  a  desert  environment,  such  as  streamflow  generated  from 
isolated,  large-magnitude  storm  events.  The  relatively  long  inter-storm  periods  effectively 
lowers  the  mean  value. 
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Figure  6.2  Monthly  temperature  time  series  at  three  elevations  and  their  spectra  before  (thin 
line)  and  after  (thick  line)  seasonal  cycle  is  removed  (cpm  is  cycles  per  month). 
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Figure  6.3  Monthly  precipitation  time  series  at  three  elevations  and  their  spectra  before 
(thin  line)  and  after  (thick  line)  seasonal  cycle  is  removed  (cpm  is  cycles  per  month). 
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Figure  6.5  Elevation  dependence  of  the  precipitation  model  parameters  in  (62). 
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6.4  P-T  Phase  Portrait  and  Animation  from  Digital  Topography 

The  phase  portrait,  or  temperature-precipitation  trajectory,  is  a  useful  tool  for 
interpretation  of  seasonal  climate  conditions.  Figure  6.6  illustrates  the  P(z,t)-T(z,t) 
trajectory  at  two  elevations  on  the  Wasatch  Front.  Plotted  on  the  horizontal  axis  is  the  first 
two  terms  of  equation  (61),  that  is  the  mean  plus  seasonal  harmonics  of  temperature.  The 
vertical  axis  represents  the  first  two  terms  of  equation  (62),  or  the  mean  plus  seasonal 
harmonics  of  precipitation.  The  dashed  line  represents  the  trajectory  of  tlie  mean 
precipitation-temperature  relation  with  elevation.  At  an  elevation  of  1,289m  (Salt  Lake 
City)  the  monthly  mean  temperature  is  10.8®C  and  precipitation  35.6  mm/mo  Note  the 
large  annual  range  in  temperature  and  the  small  range  in  precipitation  at  this  low  elevation, 
as  shown  earlier  in  Figures  4  and  5.  The  wettest  month  is  April,  likely  a  consequence  of 
the  large  scale  low  pressure  system,  the  Utah-Nevada  Low,  also  mentioned  earlier.  The 
low  pressure  system  is  reestablished  in  October,  causing  a  second  peak  of  precipitation. 
The  strong  winter  frontal  storms  certainly  have  a  very  small  contribution  at  this  low 
elevation. 

At  3,048  m  (near  the  peaks  of  the  front),  the  mean  monthly  temperature  cools  to 
0.2OC  with  mean  precipitation  increasing  to  107.  mm/mo.  Note  three  features  as  compared 
to  the  trajectory  at  1,289  m.  First,  it  has  a  smaller  temperature  and  larger  precipitation 
range.  Second,  the  direction  of  the  trajectory  is  reversed.  While  spring  is  wetter  than  fall 
at  1,289  m,  the  opposite  is  true  at  3,048  m.  Third,  the  wettest  month  is  February.  Each 
these  features  can  be  attributed  to  strong  winter  frontal  storms  and  orographic  lifting.  At 
high  elevation,  precipitation  has  a  larger  winter-summer  ratio,  a  shift  of  its  peak  toward  the 
winter  months  and  thus  a  wetter  fall  season.  The  moderately  wet  April  and  October  may 
again  reflect  the  effect  of  the  Utah-Nevada  Low  system. 

Finally,  Fan  and  Duffy  show  how  the  approach  can  be  used  in  visualization  of 
expected  temperature  and  precipitation  fields  using  digital  topography.  Each  field  is 
generated  by  applying  (61)  and/or  (62)  to  elevation  values  z(x,y)  for  a  given  month  (t).  The 
results  are  shown  as  3-d  maps  for  January,  April  and  August  in  Figure  6.7for  the  Wasatch 
Front  region.  From  the  spatial  fields  the  effect  of  elevation  on  monthly  precipitation  totals 
is  strongest  during  the  winter  months,  with  a  relatively  weak  orographic  effect  in  the 
summer  months.  Monthly  average  temperature  exhibits  large  elevation  differences  during 
the  summer  with  much  less  variation  in  the  winter. 

It  is  our  view  that  the  phase  portrait  trajectory  together  with  the  generated  fields 
from  the  DEM  serve  as  a  simple  way  of  visualizing  spatio-temporal  (seasonal)  climatic 
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patterns  in  mountainous  regions.  Extension  of  the  above  approach  to  include  aspect  will  be 
attempted  in  future  work. 


Expected  Mean  Monthly  Temperature  (*C) 


Figure  6.6  The  phase  por^t  for  the  trajectory  of  the  monthly  precipitation  versus 
temperature  at  two  elevations  and  the  trajectory  of  the  time  average  precipitation  with 
elevation. 
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Figure  6.7  January,  April,  and  August  temperature  and  precipitation  fields  generated  from 
the  empirical  model  (61)-(62)  and  digital  topography  of  the  Wasatch  Front  (after  Fan  and 
Duffy,  1993). 


CHAPTER  7.  PHILOSOPHY 


In  recent  years  national  and  international  efforts  have  focused  on  advancing  the 
fundamental  understanding  of  the  physical  processes  governing  the  circulation  of  water  in 
the  soil-groundwater-stream  system.  Today  it  is  clear  that  the  complexities  of  the  issue  go 
beyond  narrowly  defined  and  traditional  scientific  and  engineering  disciplines,  requiring  a 
multidisciplinary  approach.  From  our  perspective,  two  themes  have  emerged  which  have 
the  effect  of  unifying  the  seemingly  endless  categories  of  problems  concerning  the 
hydrologic  models  on  natural  landscapes.  The  first  is  that  applied  mathematics  and 
niimerical  and  symbolic  computation  are  the  centerpiece  of  environmental  problem  solving 
and  provide  a  framework  upon  which  experimental  pattern  and  predictability  may  emerge 
from  a  complex  environmental  system.  The  second  theme  has  to  do  with  the  question  of 
scale.  We  know  that  regularities  and  patterns  in  hydrologic  and  environmental  quality 
observations  can  be  found  at  all  scales  of  measure,  and  that  the  public  and  private  demand 
is  for  prediction  of  hydrologic  response  over  a  wide  range  of  scales:  from  soil  horizons  and 
fields  or  hillslopes,  to  geologic  formations  and  watersheds,  to  physiographic  regions  and 
whole  river  basins.  Current  predictive  models  typically  rely  on  laboratory  or  small-scale 
empirical  relations,  while  field-scale  predictions  are  made  outside  the  range  of  measured 
conditions.  That  is,  extrapolation  of  small  or  local  scale  ’laws'  (e.g.  Darcy's  law.  Pick's 
law,  solute  adsorption  isotherms  etc.)  are  assumed  to  be  representative  at  the  larger  scales. 
The  current  predictive  capability  has  served  many  useful  purposes,  but  greater  precision, 
reliability  and  integration  are  required  if  progress  is  to  be  made  towards  regional  and  global 
problems  of  water  resources  and  environmental  quality  management 

The  fundamental  challenge  of  environmentaiyhydrologic  science  and  engineering  is 
to  extend  our  predictive  capabilities  from  fluid  flows  at  the  laboratory  and  local  scale  to 
length  scales  more  pertinent  to  the  cycle  of  water  on  landscapes  such  as  individual 
watersheds  and  aquifers  (1-10's  km),  collections  of  watersheds  and  river  basins  (>100 
km).  This  opportunity  was  recently  recognized  by  the  National  Academy  of 
Sciences/National  Research  Council  report  "Opportunities  in  Hydrologic  Science",  where  it 
is  pointed  out  that  hydrologic  science  is  particularly  lacking  in  theoretical  and  experimental 
knowledge  across  scales.  In  our  opinion  this  requires  a  new  view  of  the  relation  between 
hydrologic  and  environmental  quality  data  and  models,  where  on  the  one  hand  models  and 
numerical  experiments  are  used  to  guide  the  collection  of  relevant  data,  and  on  the  other 
where  data  obtained  across  multiple  scales  is  integrated  into  a  relevant  information  base 
necessary  for  prediction  and  control.  This  approach  of  course  is  not  new  to  atmospheric 
science  where  an  international  network  of  observations  has  supported  the  critical  modeling 
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effort  of  weather  prediction.  Only  in  recent  years  has  the  US  observation-network  for 
hydrologic  and  environmental  quality  data  (e.g.  NAWDEX,  STORET,  LANDSAT,  digital 
elevation  data,...)  reached  a  point  where  integration  of  observations  and  models  can  begin 
for  the  terrestrial  component  of  the  hydrologic  cycle.  It  is  hoped  that  the  spirit  of  the 
research  conducted  through  the  present  research  grant  funded  by  ARO  to  Penn  State 
University,  reflects  this  consensus  of  the  scientific  community.  We  realize  that  the  scope 
of  this  research  effort  reported  here  is  comprehensive,  and  yet  incomplete  in  technical 
details.  We  can  only  reply  that  the  problem  of  rainfall-storage-runoff  is  comprehensive  in 
its  complexity  and  incomplete  in  its  understanding. 

7.1  Research  Directions 

The  goal  of  determining  the  essential  spatial  and  temporal  structure  of  hillslope, 
watershed,  and  river  basin  dynamics  using  low-dimensional  models,  recognizes  the 
extraordinary  computational  demands  of  modeling  the  rainfall-storage-runoff  processes 
where  topography,  soil,  geology  and  climate  variability,  necessitates  an  "effective" 
representation.  In  this  research  we  have  examined  the  practical  hypothesis  that  low¬ 
dimensional  models  (e.g.  a  few  state  variables)  are  predictive. 

A  range  of  questions  remain  concerning  the  low  dimensional  model  proposed 
above,  namely:  Under  what  conditions  are  spatial  averages  sufficient  to  'model'  the 
rainfall-storage-runoff  process?  Are  higher  spatial  moments  important?  Can  the  area- 
elevadon-aspect  approach  given  here  be  extended  to  the  case  of  random  or  fractal 
topography?  And  perhaps,  most  important,  how  are  the  flux-storage  relations 
parameterized  across  scales. 

Not  included  in  this  report  is  the  stochastic  analysis  of  hillslopes  using  spectral  and 
bispectral  theory  developed  by  Jin  and  Duffy  (1994)  as  part  of  this  research.  The 
bispectrum  and  nonlinear  time  series  analysis  paper  has  been  published  in  Water  Resources 
Research,  and  reprints  are  included.  The  second  paper  on  linear  spectral  theory  of 
hillslopes,  has  been  submitted  to  the  Journal  of  hydrology.  In  future  work  signal 
processing  and  the  interpretation  of  nonlinear  signal-response  will  be  a  major  focus  of  our 
effort.  In  addition  we  hope  to  extend  the  dynamical  model  to  include  surface  runoff, 
distributed  or  nonpoint  source  water  quality  response,  and  the  effect  of  multiple  soil  and 
geologic  strata. 
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Appendix  A 

PARAMETERIZATION  FOR  THE  RELATION  q(ii) 


Coefndents  for  q(T|)  =  ao  +  Si  iH-  82  11^+  33  are  summarized  as  follows. 


case 

ao 

a, 

32 

a, 

1 

-0.0122566 

0.0196948 

0.21227 

-0.116741 

2 

-0.00718991 

-  0.017625 

0.124908 

-0.0581907 

3 

-0.0218222 

0.115633 

0.165392 

-0.118292 

4 

-0.0360373 

0.349025 

-0.149575 

-  0.00133985 

5 

-0.00368428 

0.0640713 

0.403337 

-  0.690792 

6 

-0.0221507 

-  0.00867037 

0.151819 

-  0.0475724 

7 

-0.0104826 

0.246529 

-  0.21767 

0.061907 

8 

-0.0221879 

0.154545 

-0.153912 

0.0958753 

9 

0.0191508 

-  0.025979 

0.00798937 

0.0018399 

10 

0.623252 

-  0.282881 

0.0331941 

-  0.000339881 

11 

0.0163893 

-  0.607298 

3.36668 

-  1.75085 

12 

-0.000658114 

0.00196334 

0.00696221 

0.00452192 

13 

-0.238505 

0.780319 

2.51517 

-  1.88448 

14 

-0.732712 

0.445471 

-  0.0803704 

0.00521475 

15 

-6.14379 

4.19781 

-  0.837076 

0.0618115 

16 

-0.00680131 

0.129206 

0.110856 

-  0.0842858 

17 

-0.011539 

0.0467111 

0.192838 

-0.111842 

18 

-0.0135858 

-  0.0597895 

0.312056 

-0.151315 

19 

-0.0147494 

0.00596443 

0.223987 

-0.120729 

20 

0.00268886 

0.0213069 

0.184018 

-  0.0802037 
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Appendix  B 

PARAMETERIZATION  FOR  THE  RELATION  r(y,  t]) 


case 

fitted  polynomial 

1 

r  =  -0.600596+  1.54697  y-  1.02864  f 
+  0.653032  Ti  -  0.760397  y  ii  -  0.081 1913  ti2 

2 

r  =  -0.0295403  +  0.141486  y-  0.146076  f 
-  0.130993  Ti  +  0.117487  yri  +  0.148161  r\^ 

3 

r  =  -1.25038  +  3.02584  y-  1.92536 
+  2.3734 11  -  2.48641  yri  -  0.996918  ti2 

4 

r  =  0.25488  -  0.601226  y  +  0.363603  f 
-  1.25134  Ti  +  1.50386  yri  +  0.901529  ti2 

5 

r  =  -0.103359  +  0.649351  y  -  1.04637  72 
+  0.137549  Ti  -  0.2561 17  yii  +  0.266241  ii2 

6 

r  =  -0.448817  +  0.846004  y  -  0.40583  y2 

-  0.0013853  Ti  -  0.041432  y  ti  +  0.145737  ti2 

7 

r  s  0.894396  -  2.2241  y  +  1.39581  y2 
-  2.57258  Ti  +  3.04769  y  ii  +  1.48404  ti2 

8 

r  =  2.06444  -  5.1918  y  +  3.24758  y2 

-  3.02469 11  +  3.80659  yii  +  1.18171  ti2 

9 

r=  -27.6178  +  13.6557  y-  1.68858  f 
+  7.89528 11  -  1.95615  y  11  -  0.547036  ii2 

10 

r  =  -44.9498  +  14.2581  y-  1.13741  f 
+  8.15466 11  -  1.3029  y  11  -  0.361 1 1 1  ii2 

11 

r  =  0.856927  -  4.68138  y+  5.63224  f 
-  3.3265 11  +  10.2741  yn  +  2.94938  ti2 

12 

r  =  -0.00790433  +  0.0241201  y-  0.0198697  f- 
+  0.0037182  n  -  0.00287387  y  11  +  0.00455584  ii^ 

13 

r  =  8.47777  -  25.5956  y+  19.9379  yZ 

-  19.1595 11  +  27.2253  yii  +  10.9892  tiZ 

r  =  156.052  -  51.557  7+  4.28883  f 


-  29.8699  T)  +  4.87305  yti  +  1.43971 


14 


case 

fitted  polynomial 

15 

r  =  264.71  -  109.193  Y+  11.6957 

-  59.3592  T\  +  1 1.5059  yr\  +  3.44246  ti2 

16 

r  =  0.0175777  -  0.0178342  y-  0.0162536  y2 
-  i  .255801  Ti  +  0.440949  yt)  +  0.280376  ti2 

17 

r  =  -0.561619+  1.45122  y- 0.966134  y2 
+  0.66263  Ti  -  0.737519  yri  -  0.103081  ti2 

18 

r  = -2.18251  +  5.31406  y- 3.29844  y2 
+  2.95242  ^  -  3.52763  y  ti  -  0.921721  ti2 

19 

r  =  -J  . 32396  +  3.26534  y-  2.06584 

+  1.82922  ^  -  2.13889  yri  -  0.562959  ti2 

20 

r  =  -2.18899  +  5.04628  y-  2.90585  y2 
+  2.99582  Ti  -  3.37525  y  T]  -  0.915284  ti2 

Appendix  C 

PARAMETERIZATION  FOR  THE  RELATION  b(T|) 


Coefficients  for  b(T|)  =  ao  +  aj  t]  +  a2  11^+  di^  are  summarized  as  follows. 


case 

ao 

ai 
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.  % 

1 

-0.336185 

1.79049 

-1.89537 

0.962291 

2 

-1.10185 

4.64517 

-  5.49403 

2.39557 

3 

-0.180528 

1.1619 

-  1.07807 

0.64611 

4 

-0.0698513 

0.575338 

-0.124827 

0.215397 

5 

-0.152382 

2.95064 

-  6.23662 

9.01029 

6 

-0.50495 

1.47304 

-  1.04505 

0.318246 

7 

0.0502717 

-  0.51244 

1.02533 

-0.0774156 

8 

0.0865337 

-0.851851 

2.18505 

-  0.798875 

9 

-0.377101 

0.358386 

-  0.0760458 

0.00738185 

10 

-4.11953 

1.60227 

-  0.204627 

0.00920631 

11 

-0.233898 

1.3704 

-  1.54107 

1.01157 

12 

-0.310726 

1.68137 

-  1.57922 

0.823962 

13 

-0.338929 

1.94723 

-  1.86052 

0.97213 

16 

-0.0757071 

1.10363 

-  1.72519 

1.12782 

17 

-0.26482 

1.80463 

-  2.25957 

1.21578 

18 

-0.325671 

0.730326 

0.068719 

0.0900953 

19 

-0.4683 

2.07984 

-  2.09865 

1.02825 

20 

.-0.0111592 

0.115312 

-  0.0759964 

0.40259 

